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Introduction 

_ * 

On May 7, 1964, at the Invitation of Dr. Raymond L. Blsplinghof f , 
Associate Administrator for Advanced Research and Technology, a confer- 
ence was held in NASA Headquarters, on the subject "Heat Transfer Computer 
Programs." The conference was one of several held In Headquarters during 
the first half of 1964 for review and planning of NASA's research program 
on space vehicle thermal radiation and temperature control. 

The conference brought about a useful interchange of information 
between the NASA Centers, including JPL, on computer programs which had 
been developed for the thermal design of spacecraft. The conferees concluded 
that in addition to our annual conference in this area, a useful con- 
tinuing exchange of information could be through the establishment of a 
quarterly newsletter issued from Head luarters describing, current efforts. 

Such a newsletter, on heat transfer computer urograms, has been established 
an.l it is hoped that the publication of the present proceedings will be a 
further milestone in advancing the technology of using electronic computers 
as a tool in the thermal design of space vehicles. 

The present volume includes al*. of the formal presentations made by 
staff members of the various NAd* Centers. Their cooperation, which led 
to a highly successful conference, is ?re r tly rporecirte l. 

Coor* ’ 7, f'oolr 

T T eal usvters (Code ’7-1) 
Washington, D.C. 
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Compute r Programs Utilized in the 
Thermal Design and Analysis of Satellites 

William C. Snoddy 
MSFC 


INTRODUCTION 

Three IBM 7094 computer programs used in the thermal design of space- 
craft will be discussed. The first program is one used for integrating 
spectral reflectance or _em.it tan££. data to obtain total reflectance or 
emittance. The second program calculates the percent time a satellite 
spends in the earth r s shadow during a revolution about the earth. And the 
third program is called our General Space Thermal Program and is used to 
calculate satellite temperatures throughout a revolution about the earth. 


SPECTRAL INTEGRATION PROGRAM 

This progra»-4nt«gratre'S“ spectral radiometric data with respect to 
b lackbo< Su radiation, at any temperature or with respect to the Johnson 
solar radiation curve. One of the advantages of this program over similar 
programs is the simplicity of data input. Only those data points are 
needed which together with linear interpolation between points will define 
the spectral data curve. For example, the curves shown in figure 1 were 
entered into the computer by reading only those points marked on the 
curves. These points were chosen because it can be seen that if they are 
connected by straight lines (which the computer does when it interpolates 
linearly) the result will closely follow the original curves. Thus, for 
much spectral data it is necessary to enter only four or five points in- 
to the computer. 

Figure 2 is the computer printout for the data taken from the upper 
curve in figure 1. First is given the identification of the sample and 
then the input data is printed out. an y two of the fol lowina. (1) 

jemittance, (2) ref lec tance, . .and. (3) transmittance , the computer will cal- 
culatethe third using Kirchoff 's Law. In this case the sample was opaque 
so that the transmittance was assumed to be zero and thus the emittance 
was calculated. The integration of this data WlLlTTespect to solar-type 
radiation was of interest and therefore this integration was performed 
with respect to a blackbody curve at 5800°K and 6000°K (for other temper- 
atures it is only necessary to input the temperatures themselves) and 
also with respect to the Johnson solar curve which is stored in the 
computer.. The column on the right gives the percent of the energy outside 
the region covered by the data (in this case less than 0.32 microns and 
greater than 3.0 microns). The center three columns give the integrated 
emittance (which is equal to the absorptance by Kirchoff* 8 Law) for each 
of the three cases of interest and for various means of extrapolation to 
cover the wavelengths outside the data region. In the first of these 



columns the Integrated value is shown where a "straight" extrapolation is 
used. That is the radiometric values for the two end data points are used 
as the extrapolated values. In the next column results are given for an 
extrapolation made by fitting a curve to the data. And in the next column 
the results are shown for the case where no extrapolation is made or it can 
be thought of as an extrapolation made such that it does not alter the 
value of the integration. It might be noted that in this case, depending 
on the solar model (blackbody or Johnson) and manner of extrapolation, a 
solar absorptance ranging from 0.32 to 0.39 was obtained. 

The 7094 computer together with a S.C. 4020 also plots the data as 
shown in figure 3 (reflectance; and figure 4 (emittance or absorptance). 

The lower curve in figure 1 was integrated to obtain total emittance 
with respect to various blackbody temperatures and by mistake the Johnson 
solar curve. These results are given in figure 5 where the results for the 
Johnson curve indicates problems caused by extrapolation over lar^e regions. 
Plots of the spectral reflectance and emittance are given in figures 6 and 
7 respectively. By using the results of this program and fitting a curve 
to the total emittance as a function of temperature the variation of 
emittance with respect to temperature could easily be included in thermal 
computer programs. 


SHADOW-SUNLIGHT PR OGRAM 

In this program information regarding the passage of a satellite 
through the earth's shadow (figure 8) is obtained for use in satellite 
thermal design and analysis work. The effect that variations in the per- 
cent of the time that a satellite spends in the earth's shadow per revo- 
lution can have on temperature of a satellite is shown in figure 9 where 
a close correlation between this parameter and the internal temperature 
of Explorer VII during 320 days of orbiting can be seen. 


Great pains have been taken with this program to simplify data input 
as much as possible. For prelaunch studies it is usually convenient to 
input the data shown in figure 1C and for postlaunch calculations the 
alternate form of data input given in figure 11 may he used, since these 
parameters are usually readily available for any orbiting satellite. 

From this input the program calculates, among other things, the orbital 
positions where the satellite enters and TeaVflB the safeties sliudu m- 
and the nercant jjSaL fha «h«dr>«- Calculations are al so made for 

Successive days by extrapolation of the orbital elements consi 

p erturpatlons cause** h y &hJUfcanaaa»af -tha^jMU&i. To avoid having to 
load ~tTTe posit lorT^t he sun into the computer from the ephemerls, equations 
ware derived and Included in the program which calculate the right ascension 
and declination of the sun considering such parameters as leap years. 


Some of the results of the program are shown in figure 12. The top 
curve is the percent of time spent in the sunlight per revolution as a 
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function of days after launch for the Explorer VII satellite. This curve 
was obtained in one computer run consisting of 365 complete calculations 
and required about one minute of machine time giving a running time of 
about 1/6 second per calculation. The orbital parameters ’'ere extrapolated 
for the entire year and it was later determined that these extrapolated 
results were in good agreement with results obtained using actual orbital 
parameters. In the lower part of figure 12 is plotted additional computer 
output, namely, the orbital position of ingress into the earth's shadow, 
egress from the shadow, and perigee relative to the ascending node and as 
a function of days after launch. The position is given in terms of time 
(in minutes) required for the satellite to travel from the ascending node 
to the point of interest and an error analysis indicates accuracy to less 
than 10 seconds. This accuracy was further verified by visual observations 
of Echo I. 

This type of graph is especially useful in determining what part of 
the orbit a given tracking station will be able to "see" and thus aid in 
planning tracking and performing data analysis. This is accomplished by 
the fact that a tracking station at some fixed geographical latitude can 
at best "see" only a certain part of the satellite's orbit on any given 
day. In the case of Explorer VII and the Huntsville station, that part 
of the orbit approximately between the two .Lines on the graph in figure 12 
could be tracked at some time during the day. Thus initially Huntsville 
could track only during the daylight portions of the orbit. On about 20 
days after launch we could catch the satellite as it went into the earth's 
shadow and around 50 days after lau >ch we could track it as it was leaving 
the shadow. Knowing this information ahead of time, tracking at the 
Huntsville station was coordinated throughout the year such that special 
effort could be made to obtain data at times of particular Interest such 
as ingress, egress, and maximum and minimum sunlight conditions. 


GENERAL SPACE THERMAL PROGRAM 

MSFC has been involved in the thermal design of several spacecraft 
including Explorers I, III, IV, VII, VIII, and XI; Pioneers III and IV; 
Saturn SA-V Payload; and the Meteoroid Measurement Satellite to be flown 
on Saturn around the end of this year. Some of these spacecraft are 
shown in figure 13. It was found that weekr and even months were often 
required to set up and debug each of the computer programs used for the 
thermal design and evaluation of each of tho spacecraft. Thus was born 
the idea of a computer program which could be utilised to minimize the 
time and effort Involved in obtaining theoretical results and at the same 
time maximise the degree of sophistication which can be efficiently used 
for such calculations. 

Th e basic general equatio n used in the General Space Thermal Program 
is given in figure 14. The is a first order fourth degree differ* 

ential equ ation repres enting the lnst ant ariiWU h&At f'Iow'Info~ an3 out** of 
a aectl uiL ollhe MpAcecra tt the"t~empftratu7e^ 'if aaaumed to bm, such 

that it c an be^'pIfeT ent ed by one "effective temperature." The computer 
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will automatically set up "n" such equations depending on how many sections 
the spacecraft must be divided in order that the resulting set of equations 
will serve sufficiently as an analytical model. The first term on the 
right in the equation represents the solar flux absorbed by the section, 
the next term is the earth albedo flux absorbed, next the earth Infrared 
radiation term, then radiation away from the surface to space, internal 
heat generation, and finally conduction and radiation exchange with all 
the other sections. The bracketed factors are treated as separate functions 
which are dependent on the specific makeup of each section. 

The manner in which the program is organized is shown in figure 15. 

The program includes many special subroutines which can be used in calcula- 
ting the separate functions (A's and Q's). Each of the subroutines are 
discussed briefly below. ’ ~ 


"D" = determines whether the satellite is in the earth's shadow at 
each calculation point in the orbit 

"PS” = determines the angle between perigee and the sun in the plane 
of the orbit 

"F,;, " * calculates the geometry factor between a flat plate or the sides 
of a cylinder and the earth as a function of altitude and orien- 
tation. 

"g" - calculates the geometry factor between a sphere and the earth 
as a function of altitude 

"h" = calculates the altitude of the satellite as a function of position 
in orbit 

MAS = determines the angle between M, a vector defining the orientation 
of a section of the spacecraft, and S, a vector toward the sun 

RAM * calculates the angle between the radius vector to the satellite 
and the vector M as a function of position in orbit 

RAS * calculates the angle between the radius vector and a vector 
toward the sun 

The F subroutine might be of particular interest and its derivation 
is as follows: The geometry for planetary thermal radiation to a cylinder 

is shown in figure 16. The geometry factors for varlou' orientation and 
altitudes were calculated using a lengthly numerical progress by Ballinger, 
et.al. (references 1 and 2) and part of their results is shown in figure 17. 
In order to most efficiently enter these results into ♦■he program with an 
overall accuracy of better than two percent, a two variable curve fit was 
made such all the Information in figure 17 is contained in one fifth degree 
polynomial the coefficients of which are each fourth degree polynomials 
(figure 18). Such an efficient means of calculating the geometry factor 
allows the luxury of being able to combine the flux calculations directly 
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into the general program itself and thus reduce data handling and overall 
time required to obtain results. 

It was further found that the albedo flux could be calculated for the 
sides of a cylinder (or a flat plate) simply by multiplying the infrared 
geometry factor (F ) by the cosine of RAS and ignoring variations in the 
"azimuth" angle (figure 19) and the fact that the flux is actually a 
very complicated function of RAS. This process gives maximum errors of 
around 5% as shown in figure 20 where the dashed curves represent the results 
obtained by this method and the solid curves represent the "true" results 
from references 1 and 2. It might be noted that during a complete revolu- 
tion about the earth the errors tend to cancel out. 

T he met hod used to integrate the set of "n" first order differential 
equation is thV“Rutta -Simpson one-third rule illustrated in figure 21. 

-'"TrPorder to recta ce running time a method of varying the integration time 
step throughout a run was incorporated. This prevents the program from 
'•blowing up" and decreasing computation times by factors of 10 or more 
in some cases. 

As an example of the simplicity of using this program the case of a 
space-fixed hollow cube in orbit was considered (figure 22). The only 
information needed to set up the program is given in figure 23 and input 
data for a run is shown in figures 24 and 25. The conduction and radiation 
coefficients between sections are entered by use of a matrix (figure 25) 
and each value is given twice (i.e. Rjj » Rji) allowing the computer to 
perform a check on data input . The results for this case are shown in 
figure 26. 

Thus by use of this program it is possible to set up and obtain fairly 
sophisticated theoretical results from a computer program within 24 hours 
in many cases. 
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SATELLITE HEAT TRANSFER PROGRAMS 


by 


Edward I. Powers 

NASA Goddard Space Flight Center 
Greenbelt, Maryland 



Computer Programs for Spacecraft Thermal vlesign 
at Goddard Space Flight Center 


The thermal design of an earth satellite requires 
an accurate knowledge of both thermal radiation inputs and 
internal temperature distribution. The use of the IBM-7094 
as an analytical tool to determine these functions is now 
widely employed. 

At Goddard Space Flight Center, two programs are 
used extensively. One determines the thermal radiation 
(i.e., direct solar radiation, earth emitted radiation, 
and albedo) to a spinning flat plate. The program is geared 
toward the analysis of spin-stabilized spacecraft, but can 
be applied to non-spinning plates either earth or space 
oriented. The second program determines the temperature 
distribution and net heat flow, provided the appropriate 
conductances, radiation view factors, and thermal properties 
are known. 

A detailed discussion of these programs is presented 


below . 





THERMAL RADIATION TO A FLAT SURFACE 
ROTATING ABOUT AN ARBITRARY AXIS 
IN AN ELLIPTICAL EARTH ORBIT: APPLICATION 
TO SPIN-STABILIZED SATELLITES 



by 

Edward I. Powers 

Goddard Space Flight Center 


SUMMARY 


The derivation of total thermal radiation incident upon 
a flat plate rotating about an arbitrary axis is presented. 
The functional relationships between position in an elliptical 
earth orbit and direct solar radiation, earth-reflected solar 
radiation (albedo), and earth-emitted radiation (earthshine) 
are included. The equations have been programmed for the 
IBM 7090 digital computer, resulting in solutions which 
relate the incident radiation to spin axis orientation and 
orbital position. Several representative orbits for a typical 
geometrical configuration were analyzed and are presented 
as examples. 
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THERMAL RADIATION TO A FLAT SURFACE 
ROTATING ABOUT AN ARBITRARY AXIS 
IN AN ELLIPTICAL EARTH ORBIT: APPLICATION 
TO SPIN-STABILIZED SATELLITES 

by 

Edward I. Powers 

Goddard Space Flight Center 


INTRODUCTION 

The satisfactory operation of an artificial satellite depends upon maintaining the payload 
temperature within prescribed limits. For example, the standard batteries employed in present- 
day spacecraft generally restrict the temperature limits to 0 ' and 40 0 C. Often experiments lo- 
cated within the satellite structure further restrict this variation. 

The temperature level of a satellite may be determined by solving the instantaneous energy 
balance 


P ' Sa s A p + q. lb + q ej ^ A, ~T< + WC p 37 

where 

P = internal power dissipation, 

s = solar constant, 

a ( « solar absorptance, 

A p = instantaneous projected area for sunlight, 
q aIb * reflected solar radiation (albedo), 
q„ = earth-emitted radiation (earthshine), 

<7 * Stefan- Boltzmann constant, 

c * infrared emittance, 
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A t = total surface area, 

T 4 = mean fourth power, surface temperature, 

WC p = heat capacity of the satellite, 

^7 = time rate of change of satellite temperature. 

If the average oroital temperature is being computed, the last term is dropped. In this case al. 
heat input terms represent integrated orbital values. 

It should be noted that the above equation, as a representation of the entire satellite, is greatly 
oversimplified. The values for e and a s generally vary over the surface and great fluctuations in 
skin temperature may exist. In practice the thermal analysis consists of the development of a 
thermal model which represents a fine mesh of interconnected isothermal nodes. The appropriate 
relationship between nodes in regard to thermal conduction and radiation interchange must be 
established. 


The energy balance for each node of a thermal model may be written 


P + Sa A 

0 5 _ P 


+ q .. + a + T C (T - T ) + a V- E F (t 4 -T 4 ) 

n H nlb n nm \ m n / , nm nm \ 1 m n / 

/ \ dT n 

= «" A .. T . 4 + K ).- 3 r • 


where 

C - conductance between nodes n andm, 

E n>i = effective emissivity between nodes n andm, 

F nm = shape factor-area product between nodes n and m. 

Since the major interest at present is to determine the satellite temperature level and not the 
gradients within, a discussion of the terms in the first equation is in order. 

For most passive controlled satellites p is relatively smaii compared with the total radiation 
input and does not have a significant effect cm the satellite mean temperature. The magnitude of 
this effect, however, depends upon the « of the surface (e.g. a surface with a low absolute * ma y 
raise the internal temperature significantly because the skin has a limited capacity for 
re-radiation). 

The remaining three heat sources, direct solar heating, earth-reflected solar heating (albedo), 
and earth-emitted radiation (earthshineX represent the significant inputs to the satellite. It is 
apparent that an adequate thermal design is predicated upon a reasonably accurate knowledge of 
these thermal radiation inputs. The major source of energy— direct solar radiation— is fortunately 
the most accurately obtained. Since the sun r s rays impinging upon the satellite are virtually 



parallel, the problem is simply one of determining the instantaneous orientatior of each external 
face with respect to the solar vector. 

The calculation of earthshine is considerably less precise, requiring fundamental assumptions 
to reduce the complexity of computation. These include the assumptions that the earth is a dif- 
fusely emitting blackbody and that the surface temperature is uniform at 450°R. This permits 
direct calculation of energy input from all visible positions on the earth. Since all locations supply 
varying inputs, an integral equation must be solved to obtain the total incident energy. 

The albedo determination involves a similar integration. The earth in this case is assumed to 
be a diffusely reflecting sphere. In addition, the source intensity is a function of the satellite loca- 
tion relative to the sunlit portion of the earth. 

For a nonrotating satellite whose orientation is "fixed in space" the calculation of the thermal 
radiation .luxes is performed at any interval during the orbit. At each orbital position a particular 
flat surface may or may not "see" the entire part of the earth's cap visible from the satellite. In 
the latter case integration for albedo and earthshine must exclude the shaded portion of the cap. 
Spin-stabilized satellites which rotate uniformly about the spin axis greatly complicate the analy- 
sis. The rotation around an arbitrary axis means in general that the heat fluxes vary, since the 
visible portion of the earth is a function of this rotation. 

The analysis presented herein is based upon the energy impinging upon a flat surface whose 
orientation is defined in vector notation (by the normal vector), and which is rotating about an 
arbitrary axis. This has a much broader application than may be realized at first. Although the 
obvious application is for spin-stabilized satellites, the results apply to any body of revolution. 
Here the interest lies in the variation of flux about the axis rather than an average value per spin. 
Thus, with the orientation of a single flat plate, the impinging fluxes on cylinders and cones are 
obtained. A sphere or a shape with a variable surface curvature along its axis requires several 
or many such plates, depending upon the precision required. It can be seen that the thermal 
radiation to the entire satellite surface may be found by simply considering a handful ci appro- 
priately oriented flat plates. 

The purpose of this paper is to present, in general form, the derivation of the governing 
equations for the radiation energy sources as stated. The equations reter specifically to a flat 
surface rotating about an axis whose orientation is arbitrary. The dependence upon orbital posi- 
tion is included; the results of the numerical integration of these equations, encompassing a suit- 
able range of applicable orbits, is presented. 


COORDINATE SYSTEMS AND VECTOR REPRESENTATION 

The primary coordinate system is fixed with respect to the orbit (Figure 1). The XY plane 
lies in the plane of the orbit with the x axis coincident with the line connecting the center of the 
earth (the focus of the ellipse) and perigee. 
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X 



LEGEND 

0- EARTH CENTER 
P- PERIGEE 


As the satellite traverses the orbit, the 
instantaneous location is defined by the angle 
a. The altitude, therefore, may be determined 
at any instant by 


A(a) 


(P+jO (1 + e) 

1 + e cos a e 


(i) 


where 

A(a) = altitude, 

p » altitude at perigee. 

R = radius of the earth, 

e = eccentricity of the orbit. 


Figure 1— Coordinate systems. 


Y 



Figure 2— Orier. lion of normal and solar vectors. 


The orientation of the unit vectors (Fig- 
ure 2) representing the normals to the flat 
plates, N, and the solar vector, S, are speci- 
fied oy the following angles: 

p = angle between the projection of N on 
the xy plane and the x axis, 

y = angle between N and the Z axis, 

v - angle between the projection of S on 
the XY plane and the x axis, 

x = angle between S and tne z axis. 

The vectorial representations in the fixed 
coordinate system are thus 


N = cos P siny i + sin (S siny j + cosy k (2) 

and 

S = cos v sinX i + sinv sinX j + cos X k . (3) 

It is convenient to introduce an instantaneous coordinate system (X' Y' V ) whose origin is fixed 
in the satellite (Figure 1). The X'Y' plane lies in the orbital plane with the x' axis coincident with 
the line from the earth's center to the satellite, z* and z have the same orientation. 


The vectors N and S in the primed system are 

N = cos (/3-a) siny i' + sin (/3 - a) sinyj' + cosyk' , 
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s 


cos (v-a) sin\ i' + sin (v-a) sin\j' * cosXk' . 


( 4 ) 

(5) 



In a similar manner the expression for the spin axis vector A may be shown to be 


A - cos (S~a) sin m i' - sin ( S - a) sin ^ j ' + cos /z k ' 


( 6 ) 


where 


8 = angle between the projection of A on the XY plane and the x axis, 


m = angle between a and the Z axis. 

PLATE ROTATING ABOUT THE SPIN AXIS 

For spin- stabilized satellites the normal 
vectors N, which represent the orientation of 
the exterior surfaces, rotate about the spin 
axis. In determining both instantaneous and 
average heat fluxes, the orientation of IN as a 
function of this rotation with respect to the 
primed coordinate system (x' Y' Z' ) must be 
known. This is accomplished by defining a 
third coordinate system X" Y" z" (Figure 3). 
The x"axis is coincident with the spin axis A. 
Y" is defined so that it lies in the plane formed 
by X" and an arbitrary fixed position of the 
normal vector N o .* 

In Figure 3, the following terms may be 
evaluated: 

cos r = cos (/?-a) sin y 


Y" 



Figure 3— Double-primed coordinate system. 


(S - a) sin/z 


+ sin (/3-a) siny sin (5 -a) sin/z + cosy cos f± , 
sin r = Vl- cos J r , 

N = cost i" + sin r cos k j " + sinr sinK k “ 


(7) 

(3) 


It can be seen that the angle between the spin axis and N 0 can be computed only in the 
first and second quadrants of the X"Y" plane. This has little effect since * varies from Oto 2 77 
(a rotation), [if the vector N o lies in the third quadrant (of the X''Y" plane) the initial position for 
rotation remains in the second quadrant.] The objective now is to relate the double-primed unit 
vectors to the primed vectors in order to determine N in the primed system. 


'Since a cohylete rotation occurs, the or.entation of N e merely indicates the starting point. 


5 



Thfc projections of j " and k " in the primed coordinate system are found in the following 
manner: The value for the unit vector k" is 


m Ii. i l 


i * N 



gl * t hj 1 + / k * 
y/p, J V h* + / * 


0) 


P sin (8 “ c) sin// cos y - cos// sin (£- i) sin > , 
h " rns (/^- a) sin y - cos (8~a) sin// cosy , 

i ' cos (3 -a) sin// sin (3- a) sin y 
- sin f - a) sin// cos f/3-a) siny 

ru huIh i I v 

. ( A * N o) x A .. mi * + nj 1 + ok 1 

| * ".) x * | + n 2 + o J (^) 

where 


m r h cosm - l sin (8 -a) sin// , 


n / cos (S~a) sin// - geos/. , 


o ~ g sin (3-a) sin// - h cos (S-a) sin// . 

N may now be determined in the primed coordinate system (X 1 Y' Z') . Byrewriting the double- 
primed unit vectors, 

i “ = pi ' + qj 1 ^ r k ' ^ 


si 1 + tj ' + uk 1 


(ii) 


vi ; + wj ' + x k ' 


J 


where 


p = cos (S - a) sin// , 
q - sin (8 - a) sin// , 
r = cos fi f 
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m 


/m 2 t n 2 f o 2 


yin 2 + n 2 + o 2 


y m 2 *• n 2 + o 2 


yf 2 + h 


2 i i 2 


y4 2 ♦ h 2 


t l 


/t! 2 ♦ h 2 + / 2 

iiy substituting for i ", j", anti t" in Equation 8 


IN - ni' i bj* < ck' 


where 


(12) 


p cos 

r 


S 

sin 7 

CO S K 

h 

v s i n t 

sin k' 

q cos 

T 

+ 

t 

s i n 7 

co<; k 

+ 

w sin t 

sin k 

r cos 

7 

t 

u 

s in r 

CO S K 


x sin r 

sin k 


whore * is the angle of rotation. IN may now be evaluated at any position during rotation through- 
out the orbit in the primed system (X' Y* V ) . 


DETERMINATION OF EARTHSHINE 

The general expression for the radiation intensity dq, from a diffusely emitting source dA e , 
impinging upon a unit area, is as follows: 

I coscj cos 77 d,*. 

dq r — . ( 13 ) 

where 

1 » source intensity, 

u> - angle between the line connecting the unit area with dA e and the normal to dA e (N. ) , 
v = angle between the line connecting dA e with the unit area and the normal to the unit area, 

D = distance between the unit area and dA . 
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! iflnit- 4 — tuith') < <i(i visible hom satellite ( 4 ‘ - 0 
wl.t:n •'.nincident with f*). 


With the application of this equation to 
a flat plate of unit area at an altitude A( « > 
the radiation intensity impinging upon the 
plate is 

nl* C Pm f: 2 ” COS... co« 7 , dA f 

I J. ,l4 > 

where 

o = Stefan -Boltzmann constant, 

T r = mean black body temperature of 
the earth's surface, 

o T 4 = the heat flux leaving the source 

dA . 

C 


rioiki (hr Kt‘o»u< ti y <«1 Figure 4 (which is similar to a sketch used by Katz*) and the vector notation 
pres* Hi* o, the [tnuH in ihe equation may be defined more specifically: 


Hi* 1 R r sinfl cos</j j* + R e sin 8 sin^k 1 , 


[ A( a ) + (l"cos^)J i' + sin 6 cos </, j * *■ R c siu^ siiKpk # , 


n 2 - [A(a) + R e (l- costf)] 2 + R 2 sin 2 


0 + tan” 1 


R. sintf 


A(a) + R e ( 1 - cos 0) 


V 


» 


rns ij 


I) N 

i.»r 


| n[A(a) * R e (l-costf)] + t» R e sin# costf>t c R e sinO siiv/J 


|/[A(a) + R # ( 1 - cos 0)] J + (R e J 8in J 0 ) 

[wH-]- 


= cos’ 


*Kau, A. J,, "Determination of Thermal Radiation Incident upon the Surfaces of an Earth Satellite in an Elliptical Orbit,” tirumman 
Aircraft Engineering t orp., Rept. XP 12.20, May I960. 
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Equation 14 may now be written in the form 


< r *" 1 [*'“> ♦ «.J ( 2n r r R . sin * 

q e " it J j COS j ' 9 + tan *‘ |_A(a) + R e (l-cos<9) 

• {a [- A (a) + R t (1-cos#)] + b R # sin 6 cos<i> + c R t sin# sin$| • 

^[A(a) + R e (1-cos#)] 2 + R e 2 sin 2 #^ S2n ^ ^ ^ ’ (15) 

The limits of integration are such that the part of the earth cap visible from the instantaneous 
location of the satellite is included. For most practical cases, at least a portion of the cap is not 
visible from a plate because of its orientation. Attempts to define the appropriate integra- 
tion limits for such cases are extremely laborious. Since the equations cannot be solved 
without the aid of a high speed computer, an alternative approach is employed.* The numerical 
integration includr the entire visible part of the earth cap as stated, but the contributions of the 
elemental area that the plate does not see are deleted if the local value of cos 77 is negative. 

Equation 15 represents the flux for the instantaneous orientation of the plate. Interest also 
lies in the determination of the flux for 0 complete rotation of the plate about the satellite spin 
axis. The average value for q E is therefore 




d* 


(16) 


where * is the rotational angle. 


DETERMINATION OF ALBEDO 

The calculation for the albedo input to an orbiting plate is similar to that for the earthshine 
but involves additional basic assumptions. The reflectance of the earth depends on what is the 
visible surface— land, water, snow, cloud cover, etc. Until precise data is available which indi- 
cates the functional relationship between the albedo and the visible surface, an approximate mean 
value of 35 percent appears satisfactory. 

The reflected radiation is assumed to be diffuse. The local earth-reflected intensity varies 
with the cosine of the angle between the solar vector and the local normal. Since the input to the 
plate is a function of the source intensity of the visible elemental areas, the reflected radiation 
can be expressed by 

*Kaiz, A. J., op. cit. 
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S alb 


IT 



R, sin# 

A(a) + R ( 1 - cos #) 



A(a) + R # 


(1 - cos# 


>] 


+ b R sin 6 cos 4> + c R 

e e 


sin# sin<£ 



COS t/l 


^A(a) + R e (1 - cos 0)] 2 + R e 2 sin 2 # J 


^ 3/2 

4 *, 2 si 


in 8 dO 


where 


s = solar constant (440 BTU/hr/ft J ) 

alb = albedo factor (0.35) 

4> = angle between the solar vector S and the normal to dA e , N e , 

Cos 0 may be obtained fum: 


( 17 ) 


cos \p = S - N e 


= cos (v - a) cos# sin\ + sin (i/-a) sin# cos0 sin\ + cos\ sin# sin<£ . (18) 

The area of the earth that contributes to the albedo input is the surface common to the sunlit 
part of the earth and the part of the cap visible from the plate. The remaining area visible from 
the plate but receiving no sunlight contributes nothing to the heat input. This is accounted for in 
the numerical integration by deleting inputs when cos 4> is negative. 

The average albedo flux received by the plate per rotation about the spin axis is 


f 2 TT 

q A d« • (19) 

J 0 

CALCULATION OF DIRECT SOLAR FLUX 

The calculation of direct sunlight involves the determination of the instantaneous projected 
area of the plate with respect to the solar vector: 


q SR = S(S-N) 


= S [a cos (i > - a) sin \ + b sin (i/-a) sin A. + c cos , 


( 20 ) 


where s is the solar constant. 
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Because of its orientation the plate may or may not be facing the sun during a rotation. Neg 
ative values of S ■ IN indicate that it is not facing the sun. 


The average flux per spin is 



(21) 


Because the orientation of the plate is fixed in space, the heat flux is constant throughout the sun- 
lit portion of the orbit. To determine whether or not the satellite is within the earth's shadow at 
any orbital position requires two simple checks.* If both of the following expressions are satisfied, 
the satellite received no direct input from the sun: 


cos 0 < 0 , 

[R e + A(a)] • I sin ■ R e 

where n is the angle between the solar vector s and x , 

cos fi = S • i * , 

cosO - cos(y-a)sin\ . 


APPLICATION 

Hypothetical geometric configurations (Figure 5) have been chosen to illustrate the use 
of the equations. The external surfaces are represented by the indicated normal vectors. 

*Katz, A. J., op. cit. 





Figure 5— Geometrical configurations. 


n 



The three sketches represent three spin axis 
orientations. 

Circular orbits of 150, 550, 1050, 2000, and 
5000 statute miles, in which the solar vector lies in 
the plane of the orbit (minimum sunlight), were con- 
sidered. Figures 6-35 indicate the orbital variation 
of eaithshine and albedo.* These neat fluxes rep- 
resent mean integrated values per rotation at the 
instantaneous orbital position. The direct solar 
flux is constant in sunlight for a specific orientation. 
Appropriate values are presented, in Table 1. 
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Table 1 


Mean Solar Heat Flux per Rotation* 





b (degrees) 

P (degrees) 

• V. 

(BTU/hr/ft 1 ) 

90 

90 

0.0 

90 

45 

99.0 

90 

0 

140.0 

90 

315 

99.0 

90 

270 

0.0 

43 

45 

311 0 

45 

0 

220.0 

45 

315 

99.0 

45 

270 

0.0 

45 

225 

0.0 

0 

0 

440.0 

o ; 

315 

311.0 

0 

270 

0.0 

0 

225 

0.0 

o 

180 

0.0 


•The tiata presented represent the instantaneous day- 
light heat flux impinging upon the rotating faces of the 
configuration in Figure 5 The sun in all cases is paral- 
lel to the X axis and \ are 90 depress. 


•The mathematical model assumed the earth cap to be comprised of 5 1 2 elemental areas. The curves are based primarily on calculations 
made at 22.^ degree increments throughout the orbit Because of this, the peaks in some of the albedo curves were estimated. It should 
be noted that, because of the configurations chosen, the orbital heat fluxes in several cases are identical except for angular 
displacement. 


(Manuscript received June 26, 196)) 
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I . Summary 


An IBM 7094, three dimensional, heat transfer 
program for satellite application is discussed. The 
program may be applied to all phases of spacecraft life, 
from launch to orbital quasi-steady state. Launch phase 
heating includes inputs by radiation from the hot fairing 
and free molecule flow heating after nose cone ejection. 
Orbital heat fluxes are accounted for by direct sunlight, 
earth-reflected sunlight, and earth-emitted input to the 
satellite. 

Internal heat exchange by conduction and radiation is 
determined, provided the app opriate conductances, radiation 
view factors, and effective emittances are known. 

A simplified analysis of an active thermal controlled 
spacecraft is incorporated into the program. One or more 
surfaces may contain shutter systems which have optical 
properties specified as a function of shutter opening. The 
program evaluates the feasibility of any configuration based 
upon maintaining an internal component within its tolerable 
temperature limits. 

The program employs a matrix inversion subroutine 
(Gauss Elimination Method) to solve the heat balance 
equations. 
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II. General Discussion 

In practice, the thermal analysis of a spacecraft 
consists of the development of the thermal model which 
represents a fine mesh of interconnected isothermal nodes. 

The appropriate relationship between nodes regarding thermal 
conduction and radiation interchange must be established. 

The general energy balance for each node of a thermal 
model may be written 

P H + Sets 4 >„ * (1) 

Cr»<-T»+) r6 h A Sh -rS+ \#Cp h d% 

where 

P Q = internal power dissipation of node n 

S - solar constant 

(v g = solar absorptance 

Ap n = instantaneous projected area for sunlight of node n 
q a ib^= absorbed earth-reflected solar radiation (albedo) 

4es^ " absorbed earth-emitted radiation (earthshinei by node n 
a" = Stefan-Boltzmann constant 

e u = infrared emittance of node n 
As„ - surface area of node n 
WCp n = heat capacity of node n 
dT/dt = time rate of change of node n temperature 
C nm = thermal conductance between nodes n and m 
■ effective emittance between nodes n and m 
F nm = shape factor area product between nodes n and m 
q an ■ free molecule flow heating to node n 

q Fn = radiation from fairing absorbed by node n 

(when fairing is present external radiation 
sources are 0) 
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The solution of the n simultaneous equations 
shown above is performed by the IBM 7094 digital computer 
with a matrix inversion routine (Gauss Elimination Method) . 
However, since the inversion applies only to linear simul- 
taneous equations, the fourth power temperature terms 
must be linearized. This is done by retaining the first 
two terms in the Taylor Expansion around T . 



where T^ n is (a) the initial value of T r for steady state 
conditions or (b) the temperature at the beginning of time 
increment At during transient cases. Substituting this 
result into (1), the general equation in matrix form 
becomes 

£a*Jxf*7-/©J 

Wh&E Bhm** - Cf\tnr\ — tThrt) ~TiM ) (hfr 7)) 

fitm Vh 'f ( 3) 

D^SoisAr^n. i- fa n n, 

+ 3*6, *>„-%,++ ±.3rfat^(T7 n *-- r j) 
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Since this procedure is an approximation, the 
error involved may be significant when T in and vary 
by more than a few degrees. For steady state problems 
an iterative procedure is required to approach the correct 
value of T r . Convergence is not reached until all nodal 
temperature changes are reduced to within a specified 
tolerance. Experience has shown that if initial temper- 
ature estimates are within 50°C of the converged values, 
only 2-3 iterations are required before the results 
are virtually constant. 

Transient problems employ the same matrix inversion 
routine except each calculation represents the net energy 
exchange during a specific time interval. Since the 
solution of simultaneous equations is implicit, instability 
of the temperatures is not a problem. Thus, the results of 
high conductance and (or) low heat capacity do not adversely 
affect the numorical solution as they do in most explicit 
techniques. Care, however, must be taken in choosing a 
time interval. Although instability will not occur, no guarantee 
exists that the resulting temperatures are correct. For 
instance, calculations during a normal orbit at one minute 
intervals may be satisfactory since the environmental heat 
fluxes do not vary significantly (except perhaps at sunrise 
and sunse*' . However, when considering the launch phase, 
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calculations must be on a second basis due to the rapid 
rise and attenuation of fairing and free molecule flow 
heating rates. The short time interval is required both 
to obtain accurate integrated heat fluxes as well as to 
minimize significant temperature oscillations on space- 
craft skins and externally mounted components (lightweight). 

The same argument can be applied to sunrise and sunset 
periods during the orbital phase. The two major considera- 
tions used to determine the calculating time are, therefore: 

(1) Accuracy in obtaining integrated heat flux, and 

(2) Change in nodal temperature during a time increment. 

The effect of these inputs depends upon the particular 
requirements of the problem. 

The most straightforward procedure regarding (2) is 
to permit the program to calculate the time interval based 
upon a specified permissible temperature change. This 
method, by necessity, involves an iteration for the time 
increment determination. At present, the program provides 
only for a specified calculating time. The value is based 
primarily on previous experience gained with respect to 
the particular situation. 
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III. Discussion of Specific Environments 


A. Launch Phase 

The thermal environment encountered by the 
payload during launch consists primarily of radiation 
from the hot nose cone prior to fairing ejection and 
free molecule flow heating thereafter. Because of the 
short time involved (3 to 5 minutes) the concern is 
generally for spacecraft skins and (or) externally 
mounted components. 

Each surface element (or node) is assumed to 
view primarily one portion of the fairing whose tempera- 
ture history is known. The fairing is generally divided 
isothermally into nose cap, cone, and cylindrical sections. 

Upon ejection of the fairing, the infrared radiation 

source is replaced by aerodynamic or free molecule flow 

heating, The extent of this heating is based upon the 

3 

kinetic energy available (1/2 Pv ) of the impinging air 
and the projected area of the corresponding surface normal 
to the velocity vector. It is conservatively assumed that 
all of the kinetic energy of the air is imparted to the 
impinged surface. The analysis in this case depends upon 
a knowledge of the vehicle's trajectory in addition to the 
surface projected area. 
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B. Orbital Phase 


1. Solar Input 

The solar input may be a constant value per 
surface node or a time -dependent variable. The variation 
may be due to change in orientation with respect to the 
sun as well as eclipse by the earth. The possibility of 
internodal shading or reflections must be pre-determined 
since the program employs only net solar flux to a surface. 

2. Earthshine and Albedo 

Earthshine and albedo generally represent 
between 15% and 30% of the total energy input to a near 
earth satellite. At present , constant or average values 
per orbit are used. Albedo may be combined with the solar 
input if a transient variation is deemed necessary. 

3. Internal Power 

Internal power can be read in only as a constant. 
IV. Active Control Option 

The purpose of the active control option is to 
evaluate the design feasibility of simple shutter systems. 
One control temperature with specified limits is taken 
within the spacecraft. The effective absorptance and emit- 
tance as a function of shutter opening must be known. The 
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results indicate whether the system can handle the 
imposed extremes in thermal environment . 

V. Alternate Transient Solution 

In many instances (especially when many nodes are 
involved) the implicit matrix inversion technique becomes 
inefficient from the standpoint of computer time. A simple 
finite different explicit solution has been programmed as 
an alternative. This option proves advantageous if high 
conductances and/or low heat capacity are not present to 
restrict the calculating time interval to extremely small 
values. The small time interval is necessary to prevent 
instability of the nodal temperatures. It has been shown 
with a 25 node problem that the finite difference solution 
offers no time saving over the matrix inversion. 
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Introduction 

This paper describes b-’iefly two rather specialized 
computer programs which solve for steady-state temperatures 
over a sphere. The first program deals with opaque surfaces 
and allows for skin conduction and black body internal 
radiation. This program differs from the conventional 3-D 
heat transfer programs in that it automatically subdivides 
the sphere into a specified number of nodes and computes 
for each node, surface area, conductances, radiation shape 
factors and incident solar flux. The second progam allows 
for transmission through partially or fully transparent 
external surfaces and accounts for multiple diffuse 
reflections at the inner surfaces, but does not include 
conduction. 
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The first program was written because a number of 
GSFC satellites or parts of satellites have been spherical, 
and it has been used directly or modified to calculate 
skin temperatures of ECHO-II, and Explorer 17, and of the 
magnetometers for IMP and 0G0. The second program was 
written specifically for the proposed Rebound satellite, 
a large ECHO-type balloon which was to be constructed of 
a laminate of an aluminum mesh and a transparent plastic 
film. 

Nomenclature 

A - Surface area 

A p - Projected area 

Fi j - Shape factor from node i to node j 

H - Infrared energy incident on interior 

I - Solar energy incident on interior 

K - Conductivity 

P« - Incident albedo radiation 

P ( - Incident earth emitted radiation 

Q CO nd Conductance 

R - Radius of sphere 

s - Solar constant 

S { - Solar input 

T - Absolute temperature 
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<t> - Equatorial angle measured from the i axis 

t - Thickness of shell 

X - Fraction of surface that is opaque 

a - Internal solar absorptivity 

<* - External solar absorptivity 

e - Internal I.R. emissivity 

e - External I.R. emissivity 

P - Solar reflectivity (internal) 

P * - I.R. reflectivity (internal) 

t - Solar transmissivity 

* 

t - I.R. transmissivity 

0 - Stefan-Boltzmann constant 

9 - Polar angle of spherical segment measured 

from spin axis 

0 g - Angle between spin axis and solar vector 

Y - Angle between normal to spherical segment 

and solar vector 

Subscripts 

1 - Nodes 

j - Nodes 

p - Transmitting material 

° - Opaque material 
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Opaque Sphere Program 

The sphere is divided into a number of concentric 
rings about the spin axis as shown in Figure 1. All 
rings or nodes subtend equal angles at the origin. The 
surface area of a differential element is 

dA = R 2 sin 9 d 0 d<t> (1) 

which when integrated over a node gives 

A t = 27rR 2 (cos 0! - cos 0j) (2) 

The solar energy absorbed is proportional to the 
area of the node projected in the direction of the sun. 

The projection of an elemental area is 

dA p = dA cos v (3) 

where cos y can be determined from the dot product of the 
solar vector and the normal to the elemental area. This 
expression must be integrated for three distinct zones 
(see Figure 2) , 
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a. For 0^0^ 90 - 0 g sunlight is incident on every 
elemental area of each node. For each node the 
projection is a complete ring whose area is 

A p t = 7rR 2 cos 0 S (sin 3 0j - sin 2 0 t ) 


b. For 90 - 0 S < 0 ^ 90 + 0 g the projection is no 
longer a complete ring and it becomes necessary to 
limit the integration to sunlit areas only. The 
projected area can be demonstrated to be 

A p i * R 2 [cos 9 i (sin 2 0 g - cos 2 0 t 
- cos 0j (sin 2 0 S - cos 2 0,)^ 

- 1 , 

+ sin (cos 0i /sin 0 g ) - sin (cos 0^/sin 


+ sin 2 

9 , 

cos 0 g 

cos” 1 (-cot 

*4 

cot 

e s> 

- sin 3 

9. 

cos 9 

b 

cos” (-cot 

01 

cot 

e s )J 


c. For 90 + 0 g < 0 ^ 180 the nodes do not see the 
sun and 


A 


p t 


(4a) 


0 s) 

(4b) 


0 


(4c) 



In the expressions for the limits of 0 for each of 
the zones, it was assumed that 8 g ^ 90° . For 90 < 0 g £ 180, 
the program in effect inverts the sphere by a simple trans- 
formation of variables to a^oid defining new limits of 0 
for each zone. 

The conductance between adjacent nodes, i and i + l is 


Q 


cond 


27TKt 

tan ( - 9 - ^ +x ) 
tan 


(5) 


where the angles 0 t and 0 t + i refer to the centers of the 
two nodes. This formula can be derived by integration of 
the expression for conductance between differential areas. 

For a hollow sphere, the radiation shape factor is 


simply 
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The thermal balance at any node in steady-state is 


1 ® ^ Aj P» i + s i Aj Pgi - Aj oc j Ti + ^ Ai Fj j o€ j (Ti — T j 


V 1 


+ if Q co„d t , < T * - T < ) 


(7) 


*° solve these simultaneous equations, the terras are 
linearized by a Taylor expansion and then matrix methods 
combined with iterative procedures are used to solve for 
the temperatures. 

It is planned in the near future to publish a more 
detailed description of this program. 


Explorer 17 Thermal Model 

The Explorer 17 was not a hollow sphere and hence the 
computer program had to be modified to include the internal 
components. The majority of the electronics were mounted 
on a platform located below the sphere's equator. The 
thermal model chosen was simply that of a sphere circum- 
scribed about a cylinder as shown in Figure 3. The shape 
factors were computed by hand and read in with other data. 
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Transparent Sphere Program 

The transparent sphere program solves for the steady- 
state temperature of a partially transmitting sphere. It 
is assumed that all reflections are diffuse, that conduction 
is negligible, and that there is no diffraction as the 
energy passes through the transparent part of the sphere. 

Instead of the thermal balance consisting of one 
equation per node as in Equation 7, there are now three 
equations per node. These equations are derived by the 
Poljak^ Radiosity Method and are completely general for 
any enclosure. However, the program is specifically 
written for a sphere and hence the geometric equations 
described for the opaque sphere program are included in 
this program. 


* Eckert, E.R.G. and Drake, R. M. , Jr., "Heat and Mass 
Transfer," McGraw-Hill Book Co., 


Inc., 1959, p. 405. 
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The internal energy is divided into I.R. and visible 
(solar spectrum) components. The net solar energy incident 
internally on a node is given by 

It = Si + P 4i + ^ F i(J r XP oJ + (l-X)P p j] Ij 

0 

These equations are solved simultaneously for I t by a matrix 
so lution. 

The net I.R. energy incident internally on a node is 
given by 

H t = Pi + 2. Fi j £ Xe ° j + (I-* )e p j J Tj + 2 jFi j 

j J 

[x po J+ (l-x)p p *J Hj (9) 

These equations are solved simultaneously for H in terms 

4 

of T by a matrix solution. 

Finally, the thermal balance in steady-state for each 
node is 

r 1 ' 

ApiS [_Xao i + (1-x )c? p i J + Ai^Cxao, + (l-X)a P i) P ftt + 

[xc 0 i + (l~x)e p i"| Pi - j_x Po j + (l-X)ppi + (l-X)T-lJ Ii| 

“ Ai<7Ti ^[xe 0 t + (l-x)*p 1 3 + + (l-x)*p I 3 i + 

[xpo^i + (l-x)P^i + (1 -x)t*-iJ Hi 


( 10 ) 
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These equations are linear in and hence can be 
solved by a matrix solution to yield the steady-state 
temperature distribution. 

A complete explanation of the derivation of the 
equations for the temperature of a partially transmitting 
sphere is being planned as a NASA Technical Note for 
future publication. 



GEOMETRY OF SPHERE 



Figure 1 




EXPLORER 17 THERMAL MODEL 


spin axis 



Figure 3 


E-2570 


N66 34434 


APPLICATIONS OF THE METHOD OF MONTE CARLO TO 
PROBLEMS IN THERMAL RADIATION 
by John R. Howell 
Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio 

ABSTRACT 

A summary of the work, involving the Monte Carlo method in the so- 
lution of problems in thermal radiation transfer is presented, which 
indicates general methods previously used for solving problems in which 
radiation is coupled with other modes of energy transfer. Previous work 
involving radiation in absorbing-emitting media is included. 

An example is outlined to indicate the use of the Monte Carlo method 
in the design of a space radiator. Suggestions are given for solution of 
a complex case incorporating the effects of coupled conduction, convec- 
tion and radiation, wavelength dependent and selective surfaces, noniso- 
thermal conditions, and strongly directional or nondiffuse emitting and 
reflecting surfaces. 

A discussion is given of the factors that may affect convergence, 
running time, and accuracy of the Monte Carlo solutions and of the ad- 
vantages and disadvantages of this approach for practical problems. 

Also discussed are the case of programming for complex problems and the 
probable machine time requirements of the method. 

INTRODUCTION 

In many thermal radiation problems, integral equations of a type 
that are amenable to a Monte Carlo solution appear. One example of cur- 
rent interest from a design standpoint is a space radiator system. A 
segment of such a system might be as shown in figure 1. 

A number of analyses of such a system using conventional approaches 
(refs, 1 to 4, for example) heve been performed, but these are limited 
by certain common factors. The general formulation consists of using 
"interchange" ("shape", "view", or 'bonfiguration" ) factors, which are 
derived through an auxiliary analysis to give the fraction of energy 
leaving a given area element and arriving at another. Inherent in such 
factors is the assumption of a diffuse surface, which follows the cosine 
law for emitted and reflected energy. 
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These fractions of arriving energy are then summed by an appropriate 
integration, and the net energy flux for each area element is determined. 
Although it is possible to account for such factors as wavelength de- 
pendent properties, temperature dependent properties, and coupling with 
conduction or convection in the system studied with this type of analysis, 
it becomes quite cumbersome to evaluate the multiple integral terms that 
appear. These terms may involve integration over three dimensions and 
the wavelength interval of interest. 

Recently, a number of papers have appeared in the open literature 
applying the method of Monte Carlo to problems involving thermal radia- 
tion exchange (refs. 5 to 8). Reference 5, which is concerned with ra- 
diant transfer through absorbing-emitting media, discusses the case of 
the wavelength- and temperature-dependent radiation absorption coeffi- 
cient. The advantage of this method, as pointed out by this paper, is 
the simple evaluation of an otherwise difficult-to- solve set or in tegr al 
equations. 

In reference 6, the case of transients in the temperature of an 
absorbing-emitting gas with simultaneous heat conduction is treated ex- 
tensively. In references 7 and 8 a Monte Carlo technique is used to 
solve for the steady- state, axial heat flux distribution in a rocket 
nozzle, considering the effects of flow and convection. In reference 8 
real equilibrium propellant properties were used, and temperature, pres- 
sure, and wavelength effects on the propellant radiation absorption 
coefficient are demonstrated. 

The latter three papers all involve the solution of a set of nonlinear 
integro-differential equations, generally using a Monte Carlo technique 
for evaluation of radiation integral terms. A standard method such as 
Newton-Raphson is used for evaluation of the remaining set of equations, 
which are put in finite-difference form. 

SAMPLE MONTE CARLO PROGRAM 

In essence, the Monte Carlo technique removes the evaluations of 
multiple integrals from the computer program and replaces it with a 
simple summation. As an example of this approach, we may write the en- 
ergy balance on a volume element dVg of a typical radiator fin (fig. l) 
as 

A^T + 2 J a A,a f A,e e A,a®'a-e ^ = 2 j € A,a e A,a ^Aa (!) 

A A °A 

where the first term is the energy conducted into the volume element, 
the second term is the radiant energy absorbed by the two exposed surface 
elements of dVR (assuming symmetry about the plane of the fin), and the 
term to the right of the equality is the energy radiated from the two 
exposed area elements. For an element on a radiator tube, a convection 
term must be added. 
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If the integral term can he evaluated separately, the complete 
equation can be written in finite difference form with the integral re- 
placed by a constant for each iteration. The coupled set of equations 
that results (one for each element considered) can then be solved by any 
of a number of standard techniques. A new temperature distribution re- 
sults that can be used to reevaluate the integral terms. This procedure 
is continued until the temperatures converge, that is, until the old and 
new temperatures are the same. 

The remainder of this section will be devoted to the evaluation of 
the integral terms by the Monte Carlo procedure, leaving the methods of 
solution of the entire matrix of equations and the attendant problems of 
convergence and accuracy to be settled by the knowledge, experience, and 
intuition of the interested programmer. 

Monte Carlo Procedure 

The procedure to be presented is a simple and straightforward 
"engineering" approach to the problem. It is aimed at giving the general 
methods that can be applied with no attempt at giving the more complex 
but often faster running "shortcuts" in programming with which the liter- 
ature on Monte Carlo abounds. 

The method for this type problem reduces to this: a number 

of "bundles" of energy are assumed to be emitted from each area element 
on the radiator. The energy per bundle, C^, is taken as constant for 
all bundles within the system to simplify bookkeeping, and it can also 
be specified as a program input. The rate of bundles to be emitted from 
a given element must then be given by 


N dA = 



€ A e A dX ] dA 

°b “ 


( 2 ) 


If the direction of emission of the bundle is (t},0) as shown in 
figure 2, then for a surface with directional emissivity independent 

of wavelength, the direction of an individual bundle is given by 




cos y\ sin t] dt] 


e i“-T572 

/ cos tj sin q dt| 


( 3 ) 


6 = 2 jtRg 


( 4 ) 
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where R-^ and Rg are numbers chosen at random from a set of numbers 
evenly distributed in the range 0 < R <■ 1. The derivation of equa- 
tions (3) and (4) is given in reference 5. These equations assume an 
emissivity that is symmetrical about the normal to the surface such as 
we might expect for a real etched or plated surface. For a surface 
with a constant e^, equation (3) reduces to 

cos q - ifh (5) 

If is a nonanalytic function, equation (3) may be integrated numer- 

ically to find a curve of q as a function of R^, which may then be 
curve fit to put a function, perhaps of the form 

2 

q = A + BFt^ + <JR-^ + . . . (6) 

directly into the machine. 

Directional effects might be examined for their possible use in 
increasing radiator efficiency in this manner. Directional properties 
for a grooved surface with very strong directional effects are given in 
reference 9. Such grooving could be used to minimize interchange of 
radiation between segments of the radiator by directing the radiation 
preferentially toward space. 

Returning to the problem at hand, we may specify the wavelength 
of the bundle from the equation 


R 3* 




^ A 

e A e A dA 
\nin 

~\iax 
\nin 


6 A e A dX 


(7) 


where A min and Aj^ are the limits of the wavelength range of interest. 
For ease of use, we may again fit a curve of the form of equation (6) to 
R3 as a function of A for use in the program. 

It should be noted that e_ has been assumed independent of wave- 
length and independent of direction of emission. If this were not 

so, equation (3) would have to be carried out for many values of A, and 
equation (6) would then give a family of curves with A as the parameter 
for use in the final program. 

Similarly, a family of curves with q as parameter would result 
from equation (7). 
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It is seen from the relations derived previously that tj, d, and "h 
may be found by picking three random numbers R^, Rg, and R3 and using 
them in simple relations. If the direction of emission is known, simple 
geometry determines whether the bundle will be radiated to space or will 
strike another point on the radiator surface. For example, if the ele- 
ment of bundle emission is on the fin, it can be seen from figure 3 that 
for x taken as the distance to the center of the nearest tube, the 
energy bundle must be lost to space if 

cos n > — ( 8 ) 

x 

However, if the inequality is not satisfied, the bundle will not neces- 
sarily strike a surface. Such simple tests are useful if a large frac- 
tion of bundles are lost to space, as they save many useless calcula- 
tions of final position of the bundle. 

If the bundle does strike another surface, it may be either reflected 
or absorbed. To determine this, another number R 4 is chosen and com- 
pared to the spectral absorptivity evaluated at the bundle wavelength. 

If R 4 > a^, the bundle will be reflected; otherwise, it will be absorbed. 


If it is reflected, the new direction must be chosen again by select- 
ing at random from the correct distribution of possible angles. For a 
diffusely reflecting surface, these are given by equations (4) and (5). 

For a specular surface, we need only to find the incident angle of the 
bundle by geometry, and then the angles of reflection are given by 


^refl. ” ^incident 


1 


0 refl. - 


e + it 
incident 


( 9 ) 


For a directional surface, the direction of reflection will in general 
depend on the angle of incidence. A family of curves with the angle of 
incidence tj ' must then be specified of the form 


p rj 1 , tj cos n sin n d n 

: (10) 

p v,n cos * sin 11 dT) 

Jn* 

where p^, ^ is the reflectivity in direction t| for energy incident 

from tj * • As a practical matter, such data is scarce. One source is 
reference 10 for an idealized grooved surface. 
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The reflected bundle is new followed through any successive re- 
flections until it is either absorbed at a surface or lost to space. 

Each bundle emitted from each element is followed in the same way. 

Absorptions at each area element are tallied as they occur. Mul- 
tiplying the number of bundles finally absorbed by the energy per 
bundle gives the total radiant energy absorbed at each area element, 
and we have thus completed evaluation ctf the integral term in equation (1). 

If solar radiation is to be considered, the direction of incidence 
must be specified, and sufficient bundles are fired at our radiator to 
equal the incident solar energy. The wavelength of each incident bundle 
is again chosen by equation (8) using solar emissivities and tempera- 
tures. The bundles are followed as outlined previously. 

Running time for such a program is difficult to predict. It would 
certainly depend on many factors; however, the major ones are probably 
the astuteness of the original temperature distribution guess and the number 
of energy bundles followed. For this type of program, the running time 
would be nearly proportional to the number of bundles followed. For the 
mathematically similar but more cumbersome rocket nozzle program of ref- 
erence 8, four to five iterations took less than 10 minutes time, follow- 
ing about 75,000 energy bundles on the final few. A Newton-Raphson 
technique gave good convergence of the set of equations for each iteration. 
On the IBM 7094, the four to five iterations brought convergence to with- 
in 0.1 percent for all temperatures in the distribution. The program out- 
lined here would almost certainly run faster because the multiple ab- 
sorptions in the nozzle propellant are not present, and these absorptions 
represent a major fraction of the running time. 

Advantages and Disadvantages of the Method 

The main work involved in programming this method is the initial 
evaluation of the functions describing wavelength and directional de- 
pendence. These may involve numerical integration, which needs to be 
carried out only once. The remainder of the programming then consists 
of setting up a series of "do” loops around the siiqple algebraic equa- 
tions derived above. 

Many techniques are available (ref. 11) for estimating error in 
Monte Carlo calculations, but accuracy can usually be increased by fol- 
lowing more energy bundles (decreasing the energy per bundle). Often 
it is difficult to predict or deduce the error present in a standard mul- 
tiple numerical integration. 

The advantages of the Monte Carlo method for these programs might 
be summarized as 

(1) The ability to handle complex problems in a simple and straight- 
forward manner 
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(2) The simplicity of the resulting computer program, with the 
attendent ease of checking 

(3) The ability to control to some extent and to accurately esti- 
mate the probable errors in the solutions 

There are also certain disadvantages. It is possible that for cer- 
tain problems sufficient accuracy cannot be obtained from the straight- 
forward methods outlined without excessive computer times. More sophis- 
ticated approaches would then be required, with some loss of the advan- 
tages outlined previously. These approaches are well developed in ref- 
erence 11. 

Because of the simplicity of the equations and the independence of 
each occurrence along the path of a bundle, it is sometimes difficult to 
pinpoint or even detect the presence of program errors. Judgements made 
sure usually on the basis of the final summations, which may show little 
obvious effect of error unless gross mistakes are present. Standard 
approaches generally provide a number of possible intermediate check- 
points, lending some feeling of security to the programmer. 

CONCLUSIONS 

The Monte Carlo method may offer considerable savings in program- 
ming complexity and time for thermal radiation problems in which varia- 
tions with wavelength and direction are encountered in surface radia- 
tive properties. The method appears to offer short computer running 
times even where complex variations occur. 
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STMB0U3 

A area 

Cjj energy per Monte Carlo bundle 

Planck spectral energy distribution 
®a-e factor between area elements a and e 

k thermal conductivity 

NdA. rate of bundles emitted from a given surface element 

r outside radius of tube on tube and fin radiator 

E random number from evenly distributed set in range 0 to 1 

T absolute temperature 

V volume 

x normal distance from area element on fin to centerline of nearest 

tube 

a surface absorptivity 

e surface emissivity 

X wavelength 

angles defined in fig. 2 

P^, reflectivity in angle t| for energy incident at angle tj * 

Subscripts: 

X spectrally dependent 

a,e surface area element index, absorbing or emitting element, 

respectively 

T), 6 dependent on angle t) or 6 
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Fig. 1. - Fin-tube space radiator segment. 
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Fig. 2. - Geometry of bundle emission. 
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Fig. 3. - Criterion for loss to space 0 
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ORBITING SA1ELLITE 
PREDICTION 


SURFACE TEMEERA.TURE 
m ANALYSIS 


I. Introduction 

Temperature prediction of spacecraft orbiting the moon, earth or 
other planet 8 is an essential Manned Spacecraft Center capability. 

Existing methods i was found, however, were either too simplified or 
specialize^ for K_.v requirements. A contract was, therefore, let with 
Midwest Research Institute in Kansas City, Missouri, to develop % computer 
program to determine spacecraft heat loads and/or thi n ~; skin transient 
tempera tureF~l5Tnie~'crW moon, earth or other planets. It is this 

computer “program, which has recently been completed, that I’d like to 
introduce. The com puter program is written for the I EM 709k computer 
and is in the Fortran II language. 


II. Spacecraft Thermal Balance 

First to give some insight into what are some of the major con- 
siderations in determining the thermal balance of an orbiting spacecraft; 
figure 1 shows the principal external heat loads. The q sun represents 
solar thermal radiation coming from the sun and directly impinging 
upon the vehicle. The q albedo represents that quantity of solar 
radiation from the sun that is reflected off the planet being orbited 
and then impinging upon the orbiting vehicle. Lastly, the q planet 
is that thermal radiation from the planet orbited that impinges upon 
the vehicle. 

When the vehicle is on the far dark side of the planet, however, 
it can be seen that the vehicle would not receive any q sun or q albedo; 
therefore another consideration is the sun -shade points. 

III. Major Features and Capabilities 

Before proceeding further with the fundamentals of the program, 
however, I would like next to outline some of the major features and 
capabilities of the program. They are as follows: 

1. The program has the capability of handling up to 200 
elements; each has its own thickness, initial temperature, and thermal 
properties. ' — 


2 • Hie program also has the capability to consider eight 
different caaliMS. ajw* eight dlffe^nt ^substrate materials . TfleThennal 
properties of each may be temperature dependent".' 
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3. A spinning or fixed vehicle orientation may be considered. 
The fixed orientation nay be vith respect to the sun or the planet 
being orbited. 

k. Internal heat may be considered as a function of time. 

The program has the capability to handle eight different internal 
heat tables. 


3. Constant or variable planet temperatures nay be con* 
sidered. This feature is a most important consideration for a lunar 
orbital mission. The significance will be shown later for a hypothetical , 

lunar mission. '* 

I 

IV. Assumptions • 

In order to meet the projects objectives, it was necessary to I 

introduce cer tain simplifying a ssu mption s . I feel that it is worth- ‘ 

vhile at thiT'TEEe to Introduce ~ tKe se~ assumptions for it is up to 
each individual user or potential user to evaluate the validity of 
the assumptions for each intended application. 

1. Conduction between elements, and convection between 
the vehicle and its surroundings are neglected. When MSC feels con- 
duction between elements could be a significant factor, the program 
is used only to compute the incident heat loads and then these 
outputs are loaded as input into a transient conduction type program. 

2. No radiant interchange between elements is considered. 

3« All thermal radiation Is considered to be diffused. 

4. The vehicle's absorptivity to reflected solar radiation 
Is assumed to be equal to the vehicle's absorptivity to direct solar 
radiation. 


3. Internal heat is assumed to be uniformly distributed 
over the element. 

6. The solar constant is assumed to be independent of the 
vehicle's orbital position. 

7. On the sunlit side of a variable temperature planet, 

the planet surface temperature Is assumed to vary according to Lambert's 
cosine law. The planet temperature on the dark side of a variable 
temperature planet is, however, assumed constant. 
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V . Celestial Mechanics 

So far nothing has been said about input data. A general outline 
of the type input required shall therefore be presented next. It is 
believed that tne input is quite simple and logical. 

To computethe heat loads to an orbiting vehicle, four basic 
questions must be considered. They are: 

1. What is the location of the vehicle surface element on 
the vehicle? 

2. Where is the vehicle with respect to the planet being 

orbited? 

3- What is the celestial location of the vehicle with 
respect to space? 

k. What is the sun’s location with respect to the planet 
being orbited? 

Vehicle Coordinate System 

The first question is not applicable to a spinning satellite; 
however, for an oriented vehicle, the incident heat flux can vary 
considerably from one surface position to another. For example, 
looking at figure 2, the element on the ride of the planet oriented 
vehicle shown in the upper left hand corner receives radiation from 
the planet wheieas the element on top does not. Consequently, a 
system to answer the first question (What is the location of element 
analyzed? ) is required for the thermal analysis of an oriented 
vehicle . 

Ameaningful vehicle coordinate system is set up by replacing 
the^ c^w^tex^onfigu ^tT6h'"wrth a sphe rical mathematical model shown 
in tbe-upper right band corner oFiTgure''?: The elements on the 

sphere are selected so they have the same space orientation as the 
corresponding vehicle elements. 

For a planet or moon oriented vehicle, the surface positions on 
the sphere are defined with respect to the coordinate system illustrated 
in the lower left hand corner of the figure. For the sun oriented 
case, the surface positions on the sphere are defined with respect to 
the coordinate system illustrated in the lower right hand corner. 

For a fixed orientation, the surface elements location is defined by 
* ■ and •* • as shown and are required input data for each 
element . Ibis system can be likened to defining a position on earth 
where any position on earth can similarly be defined by giving the 
proper longitude and latitude. 
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Planet Coordinate System 


Before the Impinging heat loads emitted or reflected by the planet 
can be computed, the second question (W here la the vehicle with respect 
to the planet being orbited? ) must be answered . The answer to this 
question c&fi be obtained in terms of the planet coordinate system shown 
in figure 3 , using Kepler's equations. The input required by the 
program to c o^ute^Sterhally ''thi vehicle's position with respect to the 
planet consist only of the sdmd -major axis (a), semi -minor axis (b) of 
the orbital ellipse and the true anomaly ( 4 ') at the initial time. The 
three variables are shown in the top right hand corner of the slide which 
shows the orbit rotated into the plane of the screen, and the Xp and Yp 
axes superimposed onto the orbit plane. As shown for a spinning or fixed 
orientation, the planet coordinate system is the same. 

Celestial Coordinate System 


The third question (What is the celestial location of the vehicle 
with respect to space? ) Involves defining the or' 1 t with respect to 
space. This question can be answered in ten.u» of the celestial 
coordinate system. The vehicle's orbit can be conventionally related to 
the celestial coordinate axes by three angles which are required input 
data. These angles are shown in figure 4 where ■> »• is the right ascension 
of the ascending node, >' is the argument of perlfocus, and ■- is the 
inclination of the orbital plane. 

Any reference system could have been selected for the celestial 
coordinate system; however, the basis of selection was that it be com- 
patible with standard astronomical references to minimize input data 
compilation time and effort. With this as a criteria for selection of 
a coordinate system three systems# geocentric, mod ified 

h eliocentric and ta. describe orb i$p„about 

"The ear th , a p lanet other than earth, and the moon respectively. 

For earth and planets other than earth, the two celestial 
coordinate systems used are pictured in figure 5- 

For the moon the third celestial coordinate system selected 
is illustrated by figure 6. 

Sup's Position 


The only question that remains now to be answered is: What is 

the location of the sun with respect to the planet being orbited? 
This can be answered in.„te rms of two angles, th e right as cension 
o f the sun (RA ) a nd the dec linati on X^IX^^th angles are required 
Input data to the program and* are defined by figure 7* 


Since we were careful In our selection to choose a celestial 
coordinate system that would be compatible with standard astronomical 
references, we can obtain from the Ephemerus the RA and BBC for the 
moon, earth or any other planet for any particular date. 
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VI . Hypothetical Test Case ,. 


Case I 


A hypothetical lunar mission was run using the described program. 
For thia paxti^ lunar orbital, mij^lan*. the spa cecraft was .jplane t 
oriented an dtfa&jragiable- planet su rfac e temperature feature of the 
program was used. The mission data was selected in order that the sun 
would be very near to being in the orbital plane. 

Pertinent orbit data consisted of the following: 

Altitude 10 NM to 190 NM. 

Inclination = 13°. 


Right ascension of ascending node = . 87 °. 

Argument of perifucus 270°. 

True anomaly at initial time = 0°. 

The temperature time history of differently, oriented and painted 
surface elements’Te v eal ' several interesting characteristics. 

The results shown in figure 8 , a white element ' cools initially 
even though it is almost facing directly into the sun. A similarly 
oriented black element » shown in the same slide, however, rises to a 
peak at - 60°, then gradually drops as it turns away from the sun. 

This is true since the white element reflects a considerable amount 
of solar energy, whereas the black energy will absorb the majority 
of solar radiation that impinges upon its surface. 

In the next figure ( 9 ) the temperature curves of a black ** 
and white element facing the moon are very similar since black 
and white elements absorb long wavelength radiation almost equally. 

Also of interest shown in the figure is the hump in the black element’s 
curve at about - 100° and i = 260°. This is caused by the fact that 
both elements are briefly irradiated by the sun .just before entering 
the moon's shadow and Just after leaving its shadow. However, the 
hump shows up only in the black element's curve since it absorbs the 
solar radiation more readily. 

Comparison of both black elements in figure 10 reveals at 
- 0° the element facing the moon * is almost as hot as the element 
facing away from the moon and almost looking directly at the sun. 

This demonstrates that at low lunar orbits, the planet heat can be 
as significant i- uelar heat. 
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Cage II 

Finally, to learn the effects of not considering the moon's 
surface temperature gradient, a similar run was made of Case I, 
however, this time the moon was taken as a constant temperature planet.. 
The results shown in the final figure (ll) reveal the constant moon 
temperature" curve A is a flat curve that averages out tne maximum 
and minimum peaks of the "variable moon temperature' curve 0. As you 
can see, we have a significant variation (about 100°ft) caused by neg- 
lecting the moon’s surface temperature gradient. These resuxu, confirm 
the importance of the: variable planet. UimporuLuri: method in uoaj yiing 
lunar missions . 


Hubert A. Vogt. 

Manned opncecruft. oeni.er 




Principal External Heat Loads 



ELEMENT NO 2 ELEMENT NO 2' 


t 



MATHEMATICAL MODEL jpiANET - ORIENTED) MATHEMATICAL MODEL (SUN - OtIENTED) 

Fig. 2 - Vehicle Coordinate System for a Typical Spacecraft 





Fig. 4 - Relationship between Orbit and General 
Celestial Coordinate System 




NOTE: EARTH S X C Y C AXES LIE IN ITS EQUATORIAL PLANE 

PLANET S X c Y« AXES ARE PARALLEL TO THE ECLIPTIC PLANE 


Fig. 5 - Celestial Coordinate System for Earth (Geocentric) 
and Other Planets (Modified Heliocentric) 
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SIMPLIFIED COMPUTATION OF EARTH THERMAL AND 
ALBEDO RADIATION INCIDENT ON A SPACECRAFT 
Doyle F. Swofford 


ABSTRACT 


Hie Integral used for the computation of the earth thermal and 
albedo flux on a section of a satellite is directly integrable if the 
projected area is represented by a truncated Fourier series in the 
aspect angle. A great saving in computer time results from this rep- 
resentation, Hie series for a general surface of revolution, whose 
axis coincides with the axis of rotational symmetry of the satellite, 
is obtained by dividing the surface into elemental conical frustums 
and averaging the series for a cone over the generating curve. 
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Ap( *)) 


(oj) 


SYMBOLS 

projected, or apparent, area o. a surface viewed at an 
aspect angle n 

total area of a surface 

mean solar absorptivity of a surface 

thermal emissivity of a surface 

coefficient of the kth term in the Fourier cosine series 
for Ap/A/p( Tj) 

kth coefficient of the power series in cosine t, which 
approx traat es Ap/A^,( t) ) 

heat, input to a surface per unit time 

mean albedo of the earth 

solar radiation constant 

black-body temperature of the earth 

angular coordinates of an element of surface area on the 
earth, with the origin at the satellite 

half the angle subtended at the satellite by the earth 

semi vertex angle of a cone or frustum 

aspect angle - the angle between the spin axis of a satellite 
and the line o f view 

value of i| at the point on ear+h directly beneath the 
satellite 

angle between the earth- satellite line and the earth-sun l ? ne 
symbol, for Ap//vp 


a 


btefan-Bol t. zmnnn const ant 
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Computation of the heat Input to a spacecraft surface due to earth 
thermal and albedo radiation involves evaluation of a double integral and 
extensive tabulation of the projected area of the surface. However, a 
good closed- form approximation for the integral (requiring no numerical 
integration and no tabulation of projected areas) is obtained by substi- 
tuting a truncated Fourier series for the projected area. 

The two heat input rates are given by 

’earth *V <f e V« 0 .1 e ) 

thermal 


’albedo ?A T Sr e a s cos 3 S F < a o<V> 


where F, the so-called shape factor, is given by 


n nd Q A« 

?( W> “ij J j- <*• *> \> a 

o o 1 


da d0 


Expanding the projected area in a Fourier 6eriss 


— (n) *» i 0 + i^ cosq+ig cc.t 2 t) + i^ cos bn 


By use of trigonometric identities, this series is converted to 


7^ (»l) ■ I Q + COS 1 ] 1 2 COS^T) + 1^ COS 1 *!) 

*'*T 


In « 4 o - *2 + *4 


h » h 


To * 2^ - 
l U • 3i 4 


where 
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Taking the scalar product of the unit vectors along the satellite 
axis and the view line, ve get 


cos t) = cos cos a + sin r| e sin a cos ^ 

Ap 

Making use of *^16 relation and the series for — , ve find 

Aq, 

p ( V a o> = " cos a o)|2 I o + I 2 ein 2 n e + \ I4 sin\ e J 

+ (1 - COS 2 ^)^ COS T|J| 

+ (1 - cos^a 0 )|l2^| - sin 2 n e ^ + I4 sin 2 n e ^ - - sin 2 
+ (1 - co6"\x 0 ) Ii,^| - 2 sin 2 ^ + J stn 4 n e ^ 


— sin 2 Ti e j 


+ (l - cos-c 


Thus, the problem Is reduced to finding the coefficients 1^ in 
the expansion for Ap/Aq>( tj) . For a satellite spinning about its axis 
of symmetry, the projected areas of some surfaces are given by simple 
expressions such as these: 


*T / plate | 0 


COS T] 


o*ii*£ 

* 


~ i 1 12 p 16 k 

_ = — + — cot n + — cos c tj - cos n 

P Ia * e 2 5 jt 15* 

^cylinder • \ K" l| ‘ ^ ^ ^ ^ ^ 
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^sphere “ 4 

^hemisphere = 4 + 4 008 n 


Other surfaces are much more difficult to work with. However, since 
any surface of revolution can be divided into elemental frustums of cones, 
finding the It's of a cone as a function of its semi vertex angle 6 
will make possible a simple general solution. Dividing the cone into 
infinitesimal plane elements and integrating the series for the projected 
area of a plate around the cone from 0 to 2 n, the following result is 
obtained: 


cone 


— + — cos 2 6 - 

~ COS^6 H 

i sin 5 

5* 5* 

5 n 

2 

U» — 


COS T| 


12 3 1 * 

— - T— COS 

in in 


. 4 2 

6 + — cos & cos^ 


[- ii. ♦ cos 2 6 - cos^l 

L 15 * 3 * 3 x J 


COS^TJ 


For a general body of revolution, the projected area is found by 
summing infinitesimal frustums 

A P /'W- “ . 

%ody ” ^frust. 

[ dA 


For a generating curve , y = f(x), -evolved about the X-axis, the 
semivertex angle 6 of an infinitesimal frustum is given by 


df 

dx 


tan 6 = 
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Prom this, the 1^*6 of the infinitesimal frustum can be found as 
functions of x. Then, 


dk> 


body 



dA 



frust. 


where 


dA = 2rrf( 



dx 


In polar coordinates, revolving the curve r = r( 0 ) about the 
axis 9=0 gives 


tan 6 


dr 

r cos J + — sin 
d^ 


r sin 3 


dr 

— cos 
d0 


0 


3 


and 


dA = 2 jtt^ sin 



d0 


The surfaces of revolution will, of course, be split into sections 
small enough to justify the assumption of equal temperature throughout 
a section. 

This technique for surfaces of revolution can be applied also to 
surfaces not physically circular in cross section, bur. optically so due 
^to their spin. 



For a satellite tumbling about an axis perpendicular to the symmetry axis, 
with some residual opin about the latter, the projected area is given by 


^tumbling 


-[H 


tumbling 2rt J s P* n 


djjf 


where 


COS = Sin Tj cos 


Substituting 


^spin W ^spin ( sin 006 ^ 




we get 


tumbling " ^ Io ^spin + \ (^spin sln2r l + g (^^spin sini+T 


or 


(^tumble [ Io * 2 * 5 I| *] S pi n 




tumble ° 

^ ^ tumble ~ [j2 + f I] Jspin 

tumble = \ ^4 ^ spin 
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MONTE CARLO TECHNIQUES FOR SOLVING TRANSPORT PROBLEMS 

by Morris Perlmutter 
Lewis Research Center 

INTRODUCTION 

The Monte Carlo procedure is a model sampling technique. A model is estab- 
lished, aud the behavior of sample units in this model is followed. A sufficient 
number of sample units are followed to obtain a statistical average or macroscopic 
quantities, which are the quantities of interest. 

This technique was used in crude form by Fermi in connection with the build- 
ing of the first atomic pile. Later, Von Neumann and Ulam developed and used the 
Monte Carlo procedure extensively in developing the atomic bomb. Since then this 
technique has gained considerable use in nuclear reactor problems (refs. 1 and 2), 
and we at the Lewis Research Center have been extending it to thermal radiation 
(refs. 3 and 4), rarefied gas flows, and plasma flow problems (ref. 5). This 
technique, which requires a large number of sample histories to obtain solutions 
with small variances, is receiving greater use because of the development of the 
high-speed electronic computors. 


Method of Obtaining Sample Histories 

The sample history is the behavior of a sample unit under the conditions of 
the model. The sample unit can be an actual physical quantity such as a photon 
in thermal radiation problems or a molecule in rarefied gas flow, or it can be 
some idealized nonphysical quantity such as bundle of radiation or a bundle of 
molecules that behave according to well defined laws. For purposes of illustra- 
tion let us consider rarefied gas flow through a channel. A similar problem for 
ionized gas flow through a channel with a magnetic field across the channel is 
treated in reference 5. (The procedure applied in this sample problem would 
apply to other situations such as thermal radiation as well. ) Our sample history 
will be that of a sample molecule as it passes through the channel as shown in 
figure 1. The sample molecule is equally likely tc enter at any height x 2 be- 
tween 0 and 1 for the case of a flat plate channel. To find the height at which 
the sample molecule enters the channel, a random number is picked between 0 
and 1. This random number has an equal likelihood of being anywhere between 0 
and 1 and can be generated by use of standard techniques on the high-speed com- 
putor. If a Maxwellian gas in equilibrium in the reservoir is assumed, then the 
distribution of speeds of the molecules entering the channel is given by 

f v = ZpVe-P 2 ^ (1) 

where f y is called the probability density distribution of V. Then f y dV 
is the fraction of all the molecules that enter the channel that have a speed in 
the range V to V + dV. This depends on the parameter 3 = (2RT)"-*-/^ which 
is a function of the gas temperature T in the reservoir. 
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A velocity must now be selected for the sample molecule entering the chan- 
nel from this distribution. There' are several methods of picking from this dis- 
tribution. One method consists of pickicg two random numbers R^ and R 2 and 
setting R^ equal to V. In this case every V would be equally likely. 
Therefore, it is necessary to reject some of the values of V so that the sam- 
ple distribution will agree with f y . This is done by using this value of V 
to find f v . If f y is less than Rg, then R^ = V is rejected (fig. 2) and 
a new set of random numbers is selected. If f y is greater than Rg, then 
V-^ = R^ is used. Thus after a large number of sample molecules are picked, the 
chosen velocities will satisfy the distribution given by equation (l). There 
are many different techniques for picking from a distribution, and these are 
given in the references. In a similar manner we can pick a direction from the 
appropriate distributior. for the sample molecule entering the channel. These 
are given for a Maxwellian gas by = e/2n and f ^ = 2 cos ijr sin ijr where 
t is the cone angle and 6 is the polar angle as snown in figure I (ref. 3) . 
The molecule is followed until it has a collision with another gas molecule. 

The frequency distribution for path length to collision can be given by 

(5) 

where 2 is the path length of the sample molecule and L is the average path 
length of all the molecules at that location; the path length to the collision 
is picked from this distribution. The new molecule path must then be picked, 
and the molecule is followed in this manner until it leaves the channel. 


Obtaining the Mean Properties 


The mean properties are obtained by dividing the height of the channel into 
increments of width Axg at several locations along the length of the channel 
(fig. 2). The amount of some quantity Q that is carried by a single molecule 
across Axg at that axial position is given by 


M 



00 





( 6 ) 


where M is the mass of a molecule and 
sample molecule represents. The terms 


is the number of molecules each 


S + 

y Q and 



are the sums of the 


Q 


contained in each sample molecule passing AXg in the positive or negative 
axial direction, respectively. The number of these molecules is given by S. 
If the axial flow is desired, then Q « 1 and equation (6) reduces to 


pV x = M 




S') 


(7) 
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If the mass flow in the x 2 or vertical direction is desired, then Q = V 2 /V-^ 
and equation (6) becomes 



If the density is desired, then Q = l/V^ and equation (6) reduces to 



( 8 ) 


( 9 ) 


In this manner macroscopic flow properties in the channel can be obtained by 
summing these Q properties as the molecules pass the designated cross sections. 

In a similar manner the shear stress or drag on the wall can be calculated. 
The shear stress on the wall in the x-^ direction P x x is given by 



That is, by summing the V-^ components of velocity of all the sample molecules 
hitting an elemental area Ax-j^ on the wall (fig. 1) we can find the shear 
stress or drag at the wall at that location. If the molecules sire reflected 
diffusely, they will not contribute to the shear stress. Some results obtained 
in references 5 and 6 are shown in figures 2 to 4. 

Figure 3 shows the axial flow through a flat plate channel for a collision- 
less gas. Different channel lengths over width l were used. The solid lines 

are the analytical solutions (ref. 5), and the points are the Monte Carlo re- 

sults (ref. 6). The cross-sectional measurements were made at the entrance 
x, = 0, one quarter down the channel l / 4 and the midpoint of the channel 2/2. 
The height of the channel was divided into 20 increments. The results are sym- 
metrical around l/Z and also around x 2 = l/2. 

The densities are shown in figure 4. Again the Dack half of the channel 
is not shown because of the relation 



In figure 5 the vail shear results are presented. There ie good agreement 
between the Monte Carlo and the analytical results. 

One major drawback of this technique is the large amount of computer time 
required because of the fact chat to obtain small variations in the average 
quantities a large number of trials are needed. However, various techniques 
have been developed that allow shorter running times to be obtained. Also, the 
high-speed computers are constantly being improved. This Monte Carlo procedure 
is relatively new and will probably be used a great deal more in the future for 
solving problems on widely different subjects. It has already been found very 
useful in solving problems that would be highly intractable using mere standard 
procedures. 
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DISCUSSION OF THE GENERAL ELECTRIC PROPOSAL TO 
DEVELOP A COMPUTER PROGRAM TO SELECT THERMAL COATINGS FOR 
PASSIVE TEMPERATURE CONTROL OF SATELLITES 


by 

Robert E. Kidwell, Jr. 

NASA Goddard Space Flight Center 
Greenbelt, Maryland 



DISCUSSION OF THf crfiERAL ELECTRIC PROPOSAL TO 
DEVELOP A COMPUTER PROGRAM TO SELECT THERMAL COATINGS FOR 
PASSIVE TEMPERATURE CONTROL OF SATELLITES 

by 

Robert E. Kidwell, Jr. 

NASA Goddard Space Flight Center 
Greenbelt, Maryland 


Scope of G. E. Proposal 

The General Electric Company, Missiles and Space Division, 
submitted ja proposal j to NASA if or the development of an IBM 
7094 computer program which would select the exterior surface 
coatings for passively controlling spacecraft temperatures. | 

G. E. claims that the "trial and error" procedures currently 
used can be accomplished more rationally and can therefore 
be programmed for a digital computer. In the ASME paper, 
63-HT-41, which was presented at the ASME-AIChE Heat Transfer 
Conference at Boston in August, 1963, Costello, Harper, and 
Cline, Jdescribed the procedures] that have been used at G. E. 
j to develop a coating selection program subject to the follow- 
ing restrictions^] 

1. Steady-state conditions prevail, 

2. Heat transfer occurs by radiation only, 
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3. Temperatures are optimized at only one interior 
point in the spacecraft, and 

4. Only the solar absorptance of the external coatings 
is varied to optimize temperature. The emittance 
must remain constant at initially specified values. 

General Electric proposes to develop a generalized 
program in three steps: 

1. The program would vary both solar absorptance and 
hemispherical emittance to obtain the optimum coating patterns; 

2. Temperatures would be optimized at more than one 
interior point; and 

3. The equations would be modified to account for both 
conduction and radiation heat transfer. 

In the development of the general program, the scope 
would be restricted to steady-state heat transfer. Since the 
thermal designs of most spacecraft are based primarily on 
nearly equilibrium conditions, the proposed program could 
have wide application. An obvious extension of the proposed 
program would be to account for transient temperatures. 
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Examples 

To illustrate how the present program works, two 
examples are discussed. | The first is very simple p.nd the 
second is a little more complex. 

C ase I . - Consider a simple spacecraft which can be 
represented by equations describing the heat transfer at 
two locations or nodes. One node would represent the internal 
compartment which contains the electronic modules, batteries, etc. 
and the second node would represent the skin of the spacecraft. 

It is assumed that the internal compartment is isolated ther- 
mally from the skin, as in the case of the Vanguard Satellites, 
so that the predominant mode of heat transfer is by radiation. 

The temperature of the internal compartment is a f un ction 
of: (1) a/e of the external surfaces of the skin, and (2) 
the average per orbit of the solar and planetary heat fluxes 
incident on the skin. Consider finally three orbits that 
represent typical and extreme environments. The internal 
temperature, T, in degrees absolute to the fourth power 
is simply a linear function of the a/e of the surface coatings 
as shown in Figure 1. If T d is the desired temperature, then 
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an a/e can be selected so that the maximum difference 
4 4 

between T. and T for the three orbits is a minimum, 
d 


A selection criterion can be stated formally by defining 
a parameter A. such that 

«J 


A. = (T 4 - T 4 ) 2 

J U 


where j i*epresents a particular orbit. For a jriven a/e, 
one would obtain three values of A., one for each orbit. 

j 

The optimum coating would have an a/e such that the maximum 

value of A . for a given a/e would be a minimum, as shown in 
J 

Figure 2, where A is plotted as a function of a/e. 

The procedure for applying the selection criterion in a 
way that can be programmed for a computer can be visualized 
by referring again to Figure 2. The first step would be to 
assume a trial a/e, to compute the three A^'s and to select 
the maximum A j . The next step would be to make an incremental 
change in a/e such that the new maximum A would be less 

J 

than the previous value. The process would then be continued 
until no further decrease in the maximum A^ could be obtained. 
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Case II - Consider next a cylinder in a polar orbit 
near the earth with the axis of the cylinder parallel to 
the velocity vector and with one side of the cylinder always 
facing the earth. It is assumed that there are 10 possible 
orbital conditions based on launches at 8, 10, 12, 2, and 4 
o'clock in summer or winter. Temperature is to be controlled 
at one internal point which is radiatively coupled to the 
skin. To show the advantages of a coating pattern over one 
uniform coating, the optimum a/e's are shown in Figure 3 for 
1, 2, 4, and 8 surface nodes. Included in Figure 3 is the 
maximum deviation in °F of the internal temperature from the 
desired temperature. It is noted that when 8 surface nodes 
are used the maximum temperature deviation is reduced con- 
siderably compared to the deviation for one node. In addition, 
the coating pattern for the eight node case appears to be 
rather unusual. It would have been interesting to compare 
this solution to an independent solution by conventional 
methods using the same number of nodes. 







2 NODES 8 NODES 

AT - 32° F AT - 18° F 


Case II a/e and AT for 1 , 2, 4, and 8 

surface nodes. AT - fT - Tj) 

v d'max 


FIGURE 3 
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Comments on the Proposal 

There are a number of advantages which could be derived 

i 

from a coating selection program. The program would save the 
engineering time in the coating selection phase of the thermal 
design. In addition, ^.t would save the time lost between 
the many trial solutions that are run on the computer. For 
example, in the present trial and error procedure for coating 
selection, several steps are required. First, it is necessary 
to decide what coatings are to be applied hypothetically to 
each surface. Then it is necessary to punch new input data 
cards for the heat transfer computer program and to submit 
the deck to the computer. Without priorities over other 
work for the computer (and priorities are extremely difficult 
to obtain) several hours or more are lost before the solutions 
can be analyzed to determine what adjustments should be made 
in the coating pattern for the next trial. 

^Another advantage of a selection program is that a 
computer has infinite patience^J It can be required to compare 
and select from hundreds of numbers whereas an engineer tends 
to become fatigued. It is therefore quite probable that a 
computer could arrive at a better thermal design when there 
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are seemingly unlimited combinations to be considered. 

There are Itwo reservations about the proposed 

program. First, show does one program a computer to 

! 

perform a function in which so much intuition, judgment, 
and experience is required?/ And second, jhow much machine 
time would be required to run such a program? 

In both the ASME paper and the Proposal, none of the 
essential details were given concerning the method used to 
determine how to change the coating pattern for each trial 
so that a minimum number of steps would be required to 
converge on an optimum solution. Therefore, it is difficult 
to tell what obstacles must be overcome in developing the 
general program without actually doing most of the same 
work independently. However, jit is not nearly as difficult 
to think of ways that most of the trial and error procedure 
could be programmed for a computer for any type of satellite 
thermal design problem without necessarily achieving optimum 
coatings in an efficient 

It is difficult to know how to account for tolerances 
in the coating properties, or to decide rationally whether 
an optimum coating pattern involving delicate materials 


manner. I 
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would be preferred over a less optimum coating pattern in 
which only reliable materials would be required. 

At the present time, it is convenient to break up 
problems so that the computer time for any one run is 
limited to a few minutes. Longer runs must be made at 
night. For a 50 node problem with a'cout 24 orbital 
conditions, which is typical for Goddard satellites, the 
running time for one trial is in the order of a few minutes. 
It is obvious that a coating selection program would require 
much longer times. 

Although reservations have been expressed, jit is 
believed that the possible advantages outweigh any of 
these reservations, and that the general coating selection 
program should be developed. I 


(On June 19, 1964, contract NASw-960 was awarded to G.E. 
to develop the program) 



Spacecraft Heat Transfer Analysis 



J. A. Plamondon 
Jet Propulsion Laboratory 


INTRODUCTION 

Two somewhat related but yet distinct topics are presented below. First, 
some of our thoughts, at the Jet Propulsion Laboratory, < on the position of com- 
puting as a tool for developing the thermal design of a spacecraft are riven.] 

The purpose in doing this is to stimulate discussion on the subject of the com- 
puter bv presenting one approach to how computing might fit into thermal design. 
Second, fa new general tre heat transfer program that was rocentl” accuirec from 
Arthur D. Little, - no. is briefly described. The program is called, "Vne .lethcd 
o: Zones for the Calculation of Temperature Distribution . 11 


DISCUSSION OF COMPUTING 


Over the past few years - actually from trie earl d ys r f ’ V<-> i rorul- 
sicn Laboratory 1 s involvement in space projects - those who nave been tc 

spacecraft therrr.al cent: cl, directly or indirectly , ' ave been concerned over f ' e 
nosition that the digital computer should 'old in the overall effort o’ cace- 
craft thermal design. 

There is a general awareness of the im r ortance of the digit:! cr~'f t.c' 
ar.a of analysis, in general, in a thermal design effort. At the same "im ', i •, .. s 
also very clear oh at analysis and, in rarticula? , ccmplxcated anal s's j.:\' v nr 
the di< ital computer is not orecise, and erefere, cannot be relied urrr 1 
itself to establish the final t ermal design cf • spacecraft. Solutions c* t' : ned 
from the digital computer are simply net accurate enough, or at most, o r confi- 
dence in their accuracy is not sufficient, for teem to be the sole basts of final 
design certification. Collaboration cf the experimental and tie anal ticnl is 
essentia?, to establish confidence in the suitability and reliability cf a design. 
Nevertheless, the computer, in spite of its snortcomings , does appear tc have 
great potential, rerhaps more potential than any other single tool available to 
the temperature control engineer. 

fThe usual technioue for using computers tc solve complicated heat transfer 
problems is to lump into nodea^ zones, and so forth, incremental volumes to foim a 
network ana! og or model of the heat flow paths present in the physical rrcblem. 

The computer, of course, then analyzes this computational model.* This is fine; 
but in the process of transcribing a real problem into a computational model, the 
important ouestior arises of how well this process can be carried out for very 
complicated problems, and conseouencly, how much confidence can be placed in 
computed data? 

! The difficulties encountered in computational modeling can probablv be 
considered as arising from two sources; the modeling of the heat flow paths 



4 


PAGE 2 


present in the hardware, and the assigning of numerical values to describe the 
lumps and the interactions between lumps. For geometrically simple problems, 
the lumping or formation of a heat flow model is a relatively straightforward 
process with most computer programs having analytical guide lines. However, 
when the geometry of a problem is complicated, computational modeling becomes 
difficult and simple analytical guide lines become almost meaningless. The 
lumping process becomes partly intuitive and accurate representation becomes 
heavily dependent upon engineering understanding, Judgement, and experience. 

This leads to uncertainties and a corresponding lack of confidence in computed 
results. In the area of assigning numerical values to describe the lumps and 
their interactions, the main problem is of a different nature and centers around 
obtaining reliable thermo-physical data; for example, obtaining reliable Joint 
conductance information, i roblems in this area will certainly diminish as a 
result of current and future thermo- physical research. 

In addition to those problems associated with the modeling of hardware, 
uncertainties are also introduced because of the eouations upon which we base our 
calculations; for example, it is known that the physical idealizations which are 
generally ascribed to t"e hehaiior of radiation interchange are very liberal. 

This type of uncertainty will also diminish with future research on the physical 
behavior cf heat transport. The last type of error which is encountered in com- 
puted data results from numerical inaccuracies in computation; however, this type 
of error can be controlled since it is mathematical in nature* involving numerical 
analysis. 


Faced with uncertainties of the type just described, it is legitimate to 
ouestion the possible usages of the digital computer in thermal design and, of 
course, the consequential advantages of its utilization. 

In reviewing the development scheduling of a spacecraft, it appears that 
the use of the computer should be initiated as early as possible. Certainly, it 
should be initiated no later than the completion of preliminary design when 
mission objecti'es and general spacecraft configurations are somewhat defined in 
terms of probable size, shape, power, and so forth. The purpose would be to inves- 
tigate, as early as possible, the feasibility of preliminary design concepts as 
they relate to a particular spacecraft's configuration and mission. The results 
from such early analytical investigations do not have to be very accurate or very 
detailed; all +'nat wouxd be sought is verification that preliminary design concepts 
do indeed produce desired trends. For example, suppose that preliminarv estimates 
call for a particular vpe cf temperature control hardware such as louvers; it 
certainly would be desirable to analytically verify, early in the development 
schedule, that the louvers are really needed and, if so, that they produce the 
desired degree of control. Furthermore, and perhaps Just- as important, it would 
be very desirable to knew early in design the criticality of such a device to 
parametric, thermal variations so as to know what conditions lead to marginal and 
unreliable operation. Farlv analytical investigations in which complicated space- 
craft thermal interactions are accounted for, offer the opportunity of uncovering 
critical thermal problems so as to allow a maximum of time for their solution or 
to permit alterations in configuration and/or mission objectives with the least 
amount of disturbance. It should be pointed out that such early investigations 
can only be carried o t analytically since experimental hardware is not available, 
and that only through the use of the computer can the necessary sophistication be 
brought to bear; by this is meant the multiplicity of thermal interactions. 

Further development in the methodological use of the digital computer must be 
carried out in order tnat such early investigations can be performed with reason- 
able confidence. 
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Upon confirmation of preliminary design concepts and corrections thereof, 
the simple computational models used can then be refined to conform more and more 
closely with actual spacecraft hardware in order to perform detailed thermal 
investigations. However because of the uncertainties associated with computer 
analysis, it is essential that refinements in computational modeling be carried 
out in conjunction with experimentation, ihe inclusion of the experimental with 
the analytical forms an iterative process or procedure. With somewhat refined 
computational models, parametric studies can be undertaken to determine the 
importance of various thermal elements making up the model. For example, the 
thermal influence of a particular joint thought to be important could be analyt- 
ically studied tc determine its influence. If found to be important, an experi- 
ment could be developed to determine fairly precisely the magnitude cf the 
conductance across the joint. The experimental result would then be incorporated 
into the computational model. This sort of iteration procedure can be repeated, 
using the computational model itself to determine its own shortcomings through 
parametric study, until a fairly representative model evrlves in the sense that 
the heat transfer within a piece of hardware is well understood. Utilization cf 
the computer in this n anner aids in directing and giving credence to certain test 
programs} that is, it aids in determining just what areas are reall' important to 
test. Verification of the computational model of course recuires experimental 
correlation cf its ability to predict steadv-ctate and simcle transient conditions. 
Once correlation fas been established, it would then be possible tc investigate 
all thermal events occurring during the course of a mission, which is something 
that is entirely impractical to carry cut experimentally. Another thinp that is 
impractical experimentally, but easily done analytically, is failure mode studies 
in which various possible failures can be studied to determine their effect. 

The potential advantages of better integrating the digital computer into 
the design effort lie basically in obtaining concrete design information at an 
earlier time in the over-all development schedule. Without analvsis upon which 
we can place a high degree of confidence, concrete design information must 
necessarily wait for hardware development, fabrication, and delivery before 
experimental information can be obtained. Design changes at this point can only 
be justified when catastrophic failure or unsolvable situation is expected. Thus, 
the temperature ccntrri engineer must live with his problem with little or nr 
chance cf sidestepping it through hardware redesign. What we are seeking in 
temperature control is earlier decision capability in order to better integrate 
thermal considerations into other aspects of hardware design. To do this, we must 
be able to obtain thermal information, accurately founded, at a pace consistent 
with other design activities. There appears to oe no way, other than analysis. 

Current plans are to investigate possible ways in which + v e di;i tal com- 
puter can best be used within the framework of its current uncertainties. A 
particular configuration will be analyzed in order to define the difficulties 
encountered in using the computer and to find ways around these difficulties. 

In this investigation, testing will be used as a complement to analysis in the 
manner described previously. The object for analysis will be the Arthur D. Little 
thermal scale mcdels. The ADL models have been thoroughly tested and the accuracy 
of the results is almost unouestionable. Plans call for using the computer along 

with complementary testing to predict the thermal condition of the models. For 

this effort, none of the ADL model test results will be used as inputs for 
analysis. The' will! be used onl: as a standard for comparison with analytically 
obtained predictions. In this way, it is planned that ways of sidestepping some 

cur’wit uncerteini ies in analysis will be found. This will in turn lead tc ways 

of better exploiting the potential of the computer. 
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METHOD OF ZONES 


Next, a new heat transfer program which was recently accuired from 
Arthur D. Little, Inc. and which is called, "The Method of Zones for the Calcu- 
lation of Temperature Distribution," will be described. There are two reports 
on this program, so the description will be brief. The first report covers the 
methodology of the approach used in the program. The second part is a descrip- 
tion of how to use the program covering the input and output format and a 
program listing. 

The "Method of Zones" is an improved method for obtaining solutions to 
spacecraft associated heat transfer problems involving radiation and conduction. 
The program is one of the few general type heat transfer programs specifically 
designed for spacecraft heat transfer problems. The approach used in the program 
is not that of s irmly lumping parameters to form a nodal network analogous to an 
electrical circuit as in most general heat transfer programs. Rather, the 
approach consists of breaking a problem into zones of finite area for two- 
dimensional situations and into finite volumes for three-dimensional situations. 
The approach retains the essential characteristics of boundary value techniques 
in that each zone is characterized by prescribing conditions to the boundaries 
of each zone. The boundary conditions can represent a known heat input or out- 
put rate, or an unknown heat input or output rate expressed in terms of an 
unknown boundary temperature function such as lr. The generalized heat diffu- 
sion equation is written for each zone and is integrated over the volume of the 
zone. The integration is performed by introducing the volumetric mean tempera- 
ture of the zone. The integration of the diffusion equation results in an 
instantaneous heat balance equation which involves the heat fluxes over the 
boundaries of the zone and the rate of change in the mean temperature of the 
zone; that is, we have a differential equation expressed in terms of the time 
rate of change in mean temperature and the heat fluxes over the boundaries. The 
differential equation itself is no different than one would find in the typical 
lumped parameter program. The difference results from relating the input and 
output fluxes in the differential equation to the boundary conditions ot' the zone. 
This is done by deriving approximate formulas, based upon an escumed parabolic 
temperature distribution in the zone, in which - :he fluxes at the boundaries are 
expressed in terms of the boundary's mean temperature and the zone's mean tempera- 
ture. If the boundary conditions are ignored, the program degenerates into the 
typical lumped parameter program. tiJach zone is, thus, mathematically character- 
ized by a set of equations consisting of one differential equation for the mean 
temperature of the zone and a number of boundary equations (algebraic) for the 
temperature of the boundary. The equations are solved in a manner normal to 
other general heat transfer programs in that the differential equation is solved 
implicitly; however, the boundary equations are solved explicitly. The higher 
order of approximation resulting from the parabolic assumption permits a compli- 
cated system to be subdivided into fewer parts than is necessary when conventional 
methods are used. 

In t.h' 1 on the program, three simple examples are given to illustrate 

the application arid accuracy of the method. Criteria are also given which enable 
zone sizes to be chosen properly in practical applications. The criteria are 
derived by comma rim- ca iculaticrn made by the method of zone with exact results 
obtained for slmDle cases. 
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The numerical procedure, which is a modified Gauss-Seidel procedure, 
for solving the difference equations at each time step is discussed, and 
experience by ADL with convergence of the procedure is reviewed. The program 
also Includes checking procedures which are based on the principles of reci- 
procity and conservation of energy. These provide an initial check on the 
consistency of the input through the reciprocity principle, and a running 
check on the solution of a problem through the conservation of energy principle. 
Capability is built into the program and is controllable by the user for stopping 
the computer should the conservation of energy check show that the numerical 
solution is not being carried out to desired accuracy. 

As far as the user is concerned, the input and output formats of the 
program are as simple as any the writer has seen. However, all information 
supplied to the program must be supplied in duplicate in order that the reci- 
procity principle can check for input consistency. One of the important 
advantages in the program results from the way the equations are written. 

Being written in terms of flux Instead of temperature, permits the utilization 
of Oppenheim's network method for diffuse inter-reflections or Bobco's corollary 
for specular inter- reflections involving planar surfaces. 



ORBITER HEAT FLUX COMPUTER PROGRAM 


W.A. Hagemeyer 
Jet Propulsion Laboratory 


Abstract 


A two-pnased contract is currently underway to develop l) a set 
of parametric curves of absorbed heat flu:: for specific geometries, 
and 2^ a generalized computer program whose output is absorbed heat 
flux on the surfaces of a planet orbiting spacecraft. The absorbed 
heat fluxes include tne effects of blockage by and reflections from 
adjacent spacecraft surfaces, assuming all surfaces to be diffuse. 

The final computer program will calculate fluxes on any of ten 
arbitrarily positioned plane surfaces, with any surface properties, 
planetary radiation cnaracteristics , solar intensity, and orbital 
parameters. Outputs will be heat fluxes versus time, geometrical 
configuration factors, and radiant flux interchange factors. 


This effort started roughly in May of 1963 when JPL was in the 
midst of its Voyager Studies. Thermal Control was actively involved 
in that study, seeking to provide guidelines for basic vehicle 
configuration. The transit problem was felt to be fairly well understood, 
and effort was directed at determi:J.ng the peculiarities involved in 
making a deep space probe into an orbiter. 

At this point, it might be worth while to mention some of the 
constraints to the Voyager as they apply to the thermal design. JPL 
spacecraft are typically fairly regular polygons with the electronics 
connected to the flat faces of this polygon, which then constitutes 



what Is called the spacecraft bus. This, plus the fact that 
cylinders can be represented by flat surfaces, allowed us to restrict 
the study to flat surfaces. Mission requirements of importance were 
dual planet capability (Venus or Mars), solar panel or RTG power 
source, landing capsule, and a scientific payload which would want to 
scan the planet. Obviously a solar panel power source required sun 
orientation for the solar panels and probably the spacecraft bus, 
while an RTG would allow planet orientation. Therefore, both situations 
had to be evaluated. 

Our first estimation was that the added heat fluxes near a planet 
would be critical, especially in the case of Venus. Tnis conclusion 
was reinforced by some preliminary numbers extracted from Reference 1. 

A survey of the literature (Reference 2 ) indicated that no 
published effort had been expended in the direction of accounting 
for shading of orbiter radiator surfaces by portions of the vehicle. 

It was felt that this might be a reasonable way to eliminate some 
unwanted heat fluxes. Alternatively, this approach might allow 
adding some ueat flux during periods of spacecraft quiescence. These 
heat fluxes, in order to be successfully used, must include the 
energy reflected from the adjacent surfaces. To provide tne capability 
to more readily evaluate these effects, a contract was let to generate 
a computer program which would provide the required heat flux information. 
The required flux information will consist of the following fluxes: 
solar, planet albedo, and planetary infrared, including blockage by 
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and reflections from adjacent spacecraft surfaces. Secondary information 
output will be Incident solar, planet albedo, and planetary infrared 
fluxes, view factors, and radiant flux interchange factors. 

The contract consists of two phases. The first is a parametric 
study of particular orbiter configurations. The second is to prepare 
a generalized computer program for generating the required heat 
fluxes for any configuration of flat surfaces. 

In the parametric study, the basic configuration studied consists 
of two or three adjacent sides of a rectangular box as shown in Pig. 1, 
surface 1 being the primary radiator surface, surfaces 2 and 3 being 
secondary surfaces. Fig. 1 also delineates the dimensions used to 
describe the configuration. This configuration is then considered to 
be either sun oriented or planet oriented in three particular orbits 
as shown in Fig. 2 . Specific positions considered in each orbit are 
shown in Fig. 3. For all the above conditions, eight attitudes are 
considered varying from 100km to 30,000km, a/b and c/b (see Pig. l) 
vary from 1/4 to 1 while 0 ■ 90° and.j»£» 90°, and the surface 
properties of the primary surface and the secondary surface are fixed. 

Several of the other parameters are varied at only one altitude 
at the sub- solar point in order to indicate the trends of these 
changes. 

This parametric information will be tabulated with bounding 
values plotted. It is intended that this data will provide preliminary 
design information for the up-coming planetary orbiters. 

The generalized computer program will have as inputs any arbitrary 
set of orbital parameters, planetary characteristics, solar intensity. 


* 
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number and orientation of radiator surfaces, and a specification of 
either sun or planet orientation. Using these inputs, the program 
will generate the required heat flux Information as mentioned previously 
in addition to giving the points of entry and exit from the planet 
shadow. These heat flux outputs, as well as the radiation interchange 
factors can then be Inserted directly into a temperature computation 
program to develop a complete temperature profile for the vehicle in 
orbit. 
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BASIC APPROACH 
IN 

JPL COMPUTER PROGRAMS 


Frederick N. Magee 


Jet Propulsion Laboratory 
Pasadena, California 



Numerical analysis methods for the solution of spacecraft heat 
transfer problems were not used to a significant degree at the Jet 
Propulsion Laboratory as of March, 1962. The computer program at the 
Laboratory at that time had been prepared for another application some 
years earlier and had been found to be unsatisfactory for spacecraft 
calculations. Since it was known that considerable effort in this area 
had been invested by various industrial concerns, it seemed appropriate 
to examine and compare their methods to assure that the utmost advantage 
was taken of this experience. Uppermost in mind was the fact that the 
personnel who would derive the greatest utility from this tool were 
unacquainted with computer programs. Therefore, the absolute minimum 
of familiarization time was necessary and the routine selected should 
not be so completely automatic as to make it merely a "numbers game . 

The following criteria were established as a basis for evaluation: 

1. Siraplicity--simplicity should be sacrificed only at the 
expense of actual computation time or accuracy. Even then 
trade-offs should be carefully examined. 

2. input formats should be highly system! zed. Clear, concise 
form sheets together with a code key chart to follow. 

3. Assure that data input interpretations by the key punch 
operator are minimized. 

4. Double storage of any input data will permit external checking 
<>y •«. simple sub-routine. 

Ouct ut. format should always Include a print-out of the input 


data . 



The routine format finally chosen was that employe^ by the Hughes 
Aircraft Co. essentially due to the program simplicity and the fact that 
explicit and implicit (transient and steady-state) solutions could be 
obtained from the same input data "package". 

The input package consists of: 

1. Problem identification. 

2. Time information. 

3- Initial temperatures (T). 

4. Heat capacitance (G). 

5 . Heat generated (internal) (Q). 

6. Conductive heat connection (c). 

7. Radiation heat connection (S). 

8 . Geometry data (hov nodes are i. ter connected). 

The program capacity is of the order of $00 "nodes' and to 
date problems of 250 to 300 nodes have been run without difficulty. 

Necessity forced the development of add itions , features as our 
experience increased and specific needs were recognized. 

The transient program can now monitor the temperature of any 
desired node and output the quantity of heat required to maintain it within 
a specified temperature range. At a predetermined signal temperature a 
trigger can be activated to change certain thermal characteristics (such 
as conductivity, eraissivity, absorptivity, heat capacitance) when necessary. 

These features provide considerable latitude in examining the 
u**ny complex problems inherent in the spacecraft thermal spectrum. 

I cannot, say with great certainty that the following basic question 
has been answered, Cut here is my humble attempt: 



Why computing? 

Ultimately to be able, through analysis, t^ define the 
complete temperature spectra of a spacecraft through all phases 
of its flight within the conventional limits of engineering 
accuracy. The two advantages are greatly reduced schedule 
time and cost. It is recognized that this goal is ambitious 
but should we aspire to less we would be remiss in attending 
our professional obligation. 



